更新修改

pull/3/head
陈增辉 2025-01-14 09:25:23 +08:00
parent 9c772b8e19
commit 116a7f9bb2
9 changed files with 625 additions and 419 deletions

View File

@ -14,8 +14,8 @@
#define __CUDANVCC___ // 定义CUDA函数
#define __PRFDEBUG__
#define __PRFDEBUG_PRFINF__
//#define __PRFDEBUG__
//#define __PRFDEBUG_PRFINF__
//#define __ECHOTIMEDEBUG__
#define __TBPIMAGEDEBUG__

View File

@ -821,21 +821,32 @@ void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col
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;
void* databuffer = nullptr;
if (datetype == GDT_Float32) {
databuffer = new float[datarows * datacols];
for (int i = 0; i < datarows; i++) {
for (int j = 0; j < datacols; j++) {
((float*)databuffer)[i * datacols + j] = float(data(i, j));
}
}
}
// 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, datetype, 0, 0);
poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
databuffer, datacols, datarows, datetype, 0, 0);
}
else if (datetype == GDT_Float64) {
databuffer = new double[datarows * datacols];
for (int i = 0; i < datarows; i++) {
for (int j = 0; j < datacols; j++) {
((double*)databuffer)[i * datacols + j] = double(data(i, j));
}
}
poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows,
databuffer, datacols, datarows, datetype, 0, 0);
}
else {
}
GDALFlushCache(poDstDS);
GDALClose((GDALDatasetH)poDstDS);
// GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
@ -1171,6 +1182,53 @@ RasterExtend gdalImage::getExtend()
return extend;
}
gdalImage CreategdalImageDouble(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite, bool isEnvi)
{
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 = nullptr;
if (isEnvi) {
poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
}
else {
poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
}
GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num,
GDT_Float64, NULL); // 锟斤拷锟斤拷锟斤拷
if (need_gt) {
if (!projection.isEmpty()) {
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;
}
gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num,
Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite, bool isEnvi)
{
@ -1945,7 +2003,7 @@ bool saveEigenMatrixXd2Bin(Eigen::MatrixXd data, QString dataStrPath)
Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3);
gdalImage img=CreategdalImage(dataStrPath, data.rows(), data.cols(), 1, gt, "", false, true,true);
gdalImage img= CreategdalImageDouble(dataStrPath, data.rows(), data.cols(), 1, gt, "", false, true,true);
img.saveImage(data, 0,0,1);
@ -2862,7 +2920,7 @@ void testOutClsArr(QString filename, long* amp, long rowcount, long colcount) {
}
void testOutAntPatternTrans(QString antpatternfilename, float* antPatternArr,
void testOutAntPatternTrans(QString antpatternfilename, double* antPatternArr,
double starttheta, double deltetheta,
double startphi, double deltaphi,
long thetanum, long phinum)
@ -2889,7 +2947,7 @@ void testOutAntPatternTrans(QString antpatternfilename, float* antPatternArr,
gt(1, 2) = deltetheta;
QString antpatternfilepath = getDebugDataPath(antpatternfilename);
gdalImage ds = CreategdalImage(antpatternfilepath, thetanum, phinum, 1, gt, "", true, true, true);
gdalImage ds = CreategdalImageDouble(antpatternfilepath, thetanum, phinum, 1, gt, "", true, true, true);
ds.saveImage(antPatternMatrix, 0, 0, 1);
}

View File

@ -226,7 +226,9 @@ public:
};
// 创建影像
gdalImage CreategdalImageDouble(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false, bool isEnvi = false);
gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false, bool isEnvi = false);
gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, long espgcode, GDALDataType eType=GDT_Float32, bool need_gt = true, bool overwrite = false,bool isENVI=false);
gdalImageComplex CreategdalImageComplex(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false);
@ -273,7 +275,7 @@ bool saveEigenMatrixXd2Bin(Eigen::MatrixXd data, QString dataStrPath);
// 测试
void testOutAntPatternTrans(QString antpatternfilename, float* antPatternArr, double starttheta, double deltetheta, double startphi, double deltaphi, long thetanum, long phinum);
void testOutAntPatternTrans(QString antpatternfilename, double* antPatternArr, double starttheta, double deltetheta, double startphi, double deltaphi, long thetanum, long phinum);
void testOutAmpArr(QString filename, float* amp, long rowcount, long colcount);
void testOutAmpArr(QString filename, double* amp, long rowcount, long colcount);

View File

@ -18,64 +18,64 @@
__device__ float GPU_getSigma0dB(CUDASigmaParam param, float theta) {//线性值
float sigma = param.p1 + param.p2 * exp(-param.p3 * theta) + param.p4 * cos(param.p5 * theta + param.p6);
__device__ double GPU_getSigma0dB(CUDASigmaParam param, double theta) {//线性值
double sigma = param.p1 + param.p2 * exp(-param.p3 * theta) + param.p4 * cos(param.p5 * theta + param.p6);
return sigma;
}
__device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
float RstX, float RstY, float RstZ,
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ
double RstX, double RstY, double RstZ,
double antXaxisX, double antXaxisY, double antXaxisZ,
double antYaxisX, double antYaxisY, double antYaxisZ,
double antZaxisX, double antZaxisY, double antZaxisZ,
double antDirectX, double antDirectY, double antDirectZ
) {
CUDAVectorEllipsoidal result{ 0,0,-1 };
// 求解天线增益
float Xst = -1 * RstX; // 卫星 --> 地面
float Yst = -1 * RstY;
float Zst = -1 * RstZ;
float AntXaxisX = antXaxisX;
float AntXaxisY = antXaxisY;
float AntXaxisZ = antXaxisZ;
float AntYaxisX = antYaxisX;
float AntYaxisY = antYaxisY;
float AntYaxisZ = antYaxisZ;
float AntZaxisX = antZaxisX;
float AntZaxisY = antZaxisY;
float AntZaxisZ = antZaxisZ;
double Xst = -1 * RstX; // 卫星 --> 地面
double Yst = -1 * RstY;
double Zst = -1 * RstZ;
double AntXaxisX = antXaxisX;
double AntXaxisY = antXaxisY;
double AntXaxisZ = antXaxisZ;
double AntYaxisX = antYaxisX;
double AntYaxisY = antYaxisY;
double AntYaxisZ = antYaxisZ;
double AntZaxisX = antZaxisX;
double AntZaxisY = antZaxisY;
double AntZaxisZ = antZaxisZ;
// 归一化
float RstNorm = sqrtf(Xst * Xst + Yst * Yst + Zst * Zst);
float AntXaxisNorm = sqrtf(AntXaxisX * AntXaxisX + AntXaxisY * AntXaxisY + AntXaxisZ * AntXaxisZ);
float AntYaxisNorm = sqrtf(AntYaxisX * AntYaxisX + AntYaxisY * AntYaxisY + AntYaxisZ * AntYaxisZ);
float AntZaxisNorm = sqrtf(AntZaxisX * AntZaxisX + AntZaxisY * AntZaxisY + AntZaxisZ * AntZaxisZ);
double RstNorm = sqrtf(Xst * Xst + Yst * Yst + Zst * Zst);
double AntXaxisNorm = sqrtf(AntXaxisX * AntXaxisX + AntXaxisY * AntXaxisY + AntXaxisZ * AntXaxisZ);
double AntYaxisNorm = sqrtf(AntYaxisX * AntYaxisX + AntYaxisY * AntYaxisY + AntYaxisZ * AntYaxisZ);
double AntZaxisNorm = sqrtf(AntZaxisX * AntZaxisX + AntZaxisY * AntZaxisY + AntZaxisZ * AntZaxisZ);
float Rx = Xst / RstNorm;
float Ry = Yst / RstNorm;
float Rz = Zst / RstNorm;
float Xx = AntXaxisX / AntXaxisNorm;
float Xy = AntXaxisY / AntXaxisNorm;
float Xz = AntXaxisZ / AntXaxisNorm;
float Yx = AntYaxisX / AntYaxisNorm;
float Yy = AntYaxisY / AntYaxisNorm;
float Yz = AntYaxisZ / AntYaxisNorm;
float Zx = AntZaxisX / AntZaxisNorm;
float Zy = AntZaxisY / AntZaxisNorm;
float Zz = AntZaxisZ / AntZaxisNorm;
double Rx = Xst / RstNorm;
double Ry = Yst / RstNorm;
double Rz = Zst / RstNorm;
double Xx = AntXaxisX / AntXaxisNorm;
double Xy = AntXaxisY / AntXaxisNorm;
double Xz = AntXaxisZ / AntXaxisNorm;
double Yx = AntYaxisX / AntYaxisNorm;
double Yy = AntYaxisY / AntYaxisNorm;
double Yz = AntYaxisZ / AntYaxisNorm;
double Zx = AntZaxisX / AntZaxisNorm;
double Zy = AntZaxisY / AntZaxisNorm;
double Zz = AntZaxisZ / AntZaxisNorm;
float Xant = (Rx * Yy * Zz - Rx * Yz * Zy - Ry * Yx * Zz + Ry * Yz * Zx + Rz * Yx * Zy - Rz * Yy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
float Yant = -(Rx * Xy * Zz - Rx * Xz * Zy - Ry * Xx * Zz + Ry * Xz * Zx + Rz * Xx * Zy - Rz * Xy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
float Zant = (Rx * Xy * Yz - Rx * Xz * Yy - Ry * Xx * Yz + Ry * Xz * Yx + Rz * Xx * Yy - Rz * Xy * Yx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
double Xant = (Rx * Yy * Zz - Rx * Yz * Zy - Ry * Yx * Zz + Ry * Yz * Zx + Rz * Yx * Zy - Rz * Yy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
double Yant = -(Rx * Xy * Zz - Rx * Xz * Zy - Ry * Xx * Zz + Ry * Xz * Zx + Rz * Xx * Zy - Rz * Xy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
double Zant = (Rx * Xy * Yz - Rx * Xz * Yy - Ry * Xx * Yz + Ry * Xz * Yx + Rz * Xx * Yy - Rz * Xy * Yx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
// 计算theta 与 phi
float Norm = sqrtf(Xant * Xant + Yant * Yant + Zant * Zant); // 计算 pho
float ThetaAnt = acosf(Zant / Norm); // theta 与 Z轴的夹角
float PhiAnt = atanf(Yant / Xant); // -pi/2 ~pi/2
double Norm = sqrtf(Xant * Xant + Yant * Yant + Zant * Zant); // 计算 pho
double ThetaAnt = acosf(Zant / Norm); // theta 与 Z轴的夹角
double PhiAnt = atanf(Yant / Xant); // -pi/2 ~pi/2
if (abs(Yant) < PRECISIONTOLERANCE) { // X轴上
@ -112,20 +112,20 @@ __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
return result;
}
__device__ float GPU_BillerInterpAntPattern(float* antpattern,
float starttheta, float startphi, float dtheta, float dphi,
__device__ double GPU_BillerInterpAntPattern(double* antpattern,
double starttheta, double startphi, double dtheta, double dphi,
long thetapoints, long phipoints,
float searththeta, float searchphi) {
float stheta = searththeta;
float sphi = searchphi;
double searththeta, double searchphi) {
double stheta = searththeta;
double sphi = searchphi;
if (stheta > 90) {
return 0;
}
else {}
float pthetaid = (stheta - starttheta) / dtheta;//
float pphiid = (sphi - startphi) / dphi;
double pthetaid = (stheta - starttheta) / dtheta;//
double pphiid = (sphi - startphi) / dphi;
long lasttheta = floorf(pthetaid);
long nextTheta = lasttheta + 1;
@ -139,18 +139,18 @@ __device__ float GPU_BillerInterpAntPattern(float* antpattern,
return 0;
}
else {
float x = stheta;
float y = sphi;
double x = stheta;
double y = sphi;
float x1 = lasttheta * dtheta + starttheta;
float x2 = nextTheta * dtheta + starttheta;
float y1 = lastphi * dphi + startphi;
float y2 = nextPhi * dphi + startphi;
double x1 = lasttheta * dtheta + starttheta;
double x2 = nextTheta * dtheta + starttheta;
double y1 = lastphi * dphi + startphi;
double y2 = nextPhi * dphi + startphi;
float z11 = antpattern[lasttheta * phipoints + lastphi];
float z12 = antpattern[lasttheta * phipoints + nextPhi];
float z21 = antpattern[nextTheta * phipoints + lastphi];
float z22 = antpattern[nextTheta * phipoints + nextPhi];
double z11 = antpattern[lasttheta * phipoints + lastphi];
double z12 = antpattern[lasttheta * phipoints + nextPhi];
double z21 = antpattern[nextTheta * phipoints + lastphi];
double z22 = antpattern[nextTheta * phipoints + nextPhi];
//z11 = powf(10, z11 / 10); // dB-> 线性
@ -158,7 +158,7 @@ __device__ float GPU_BillerInterpAntPattern(float* antpattern,
//z21 = powf(10, z21 / 10);
//z22 = powf(10, z22 / 10);
float GainValue = (z11 * (x2 - x) * (y2 - y)
double GainValue = (z11 * (x2 - x) * (y2 - y)
+ z21 * (x - x1) * (y2 - y)
+ z12 * (x2 - x) * (y - y1)
+ z22 * (x - x1) * (y - y1));
@ -167,12 +167,12 @@ __device__ float GPU_BillerInterpAntPattern(float* antpattern,
}
}
__device__ cuComplex GPU_calculationEcho(float sigma0, float TransAnt, float ReciveAnt,
float localangle, float R, float slopeangle, float Pt, float lamda) {
float amp = Pt * TransAnt * ReciveAnt;
__device__ cuComplex GPU_calculationEcho(double sigma0, double TransAnt, double ReciveAnt,
double localangle, double R, double slopeangle, double Pt, double lamda) {
double amp = Pt * TransAnt * ReciveAnt;
amp = amp * sigma0;
amp = amp / (powf(4 * LAMP_CUDA_PI, 2) * powf(R, 4)); // 反射强度
float phi = (-4 * LAMP_CUDA_PI / lamda) * R;
double phi = (-4 * LAMP_CUDA_PI / lamda) * R;
cuComplex echophi = make_cuComplex(0, phi);
cuComplex echophiexp = cuCexpf(echophi);
cuComplex echo = make_cuComplex(echophiexp.x * amp, echophiexp.y * amp);
@ -180,57 +180,57 @@ __device__ cuComplex GPU_calculationEcho(float sigma0, float TransAnt, float Re
}
__global__ void CUDA_SatelliteAntDirectNormal(float* RstX, float* RstY, float* RstZ,
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ,
float* thetaAnt, float* phiAnt
__global__ void CUDA_SatelliteAntDirectNormal(double* RstX, double* RstY, double* RstZ,
double antXaxisX, double antXaxisY, double antXaxisZ,
double antYaxisX, double antYaxisY, double antYaxisZ,
double antZaxisX, double antZaxisY, double antZaxisZ,
double antDirectX, double antDirectY, double antDirectZ,
double* thetaAnt, double* phiAnt
, long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float Xst = -1 * RstX[idx]; // 卫星 --> 地面
float Yst = -1 * RstY[idx];
float Zst = -1 * RstZ[idx];
float AntXaxisX = antXaxisX;
float AntXaxisY = antXaxisY;
float AntXaxisZ = antXaxisZ;
float AntYaxisX = antYaxisX;
float AntYaxisY = antYaxisY;
float AntYaxisZ = antYaxisZ;
float AntZaxisX = antZaxisX;
float AntZaxisY = antZaxisY;
float AntZaxisZ = antZaxisZ;
double Xst = -1 * RstX[idx]; // 卫星 --> 地面
double Yst = -1 * RstY[idx];
double Zst = -1 * RstZ[idx];
double AntXaxisX = antXaxisX;
double AntXaxisY = antXaxisY;
double AntXaxisZ = antXaxisZ;
double AntYaxisX = antYaxisX;
double AntYaxisY = antYaxisY;
double AntYaxisZ = antYaxisZ;
double AntZaxisX = antZaxisX;
double AntZaxisY = antZaxisY;
double AntZaxisZ = antZaxisZ;
// 归一化
float RstNorm = sqrtf(Xst * Xst + Yst * Yst + Zst * Zst);
float AntXaxisNorm = sqrtf(AntXaxisX * AntXaxisX + AntXaxisY * AntXaxisY + AntXaxisZ * AntXaxisZ);
float AntYaxisNorm = sqrtf(AntYaxisX * AntYaxisX + AntYaxisY * AntYaxisY + AntYaxisZ * AntYaxisZ);
float AntZaxisNorm = sqrtf(AntZaxisX * AntZaxisX + AntZaxisY * AntZaxisY + AntZaxisZ * AntZaxisZ);
double RstNorm = sqrtf(Xst * Xst + Yst * Yst + Zst * Zst);
double AntXaxisNorm = sqrtf(AntXaxisX * AntXaxisX + AntXaxisY * AntXaxisY + AntXaxisZ * AntXaxisZ);
double AntYaxisNorm = sqrtf(AntYaxisX * AntYaxisX + AntYaxisY * AntYaxisY + AntYaxisZ * AntYaxisZ);
double AntZaxisNorm = sqrtf(AntZaxisX * AntZaxisX + AntZaxisY * AntZaxisY + AntZaxisZ * AntZaxisZ);
float Rx = Xst / RstNorm;
float Ry = Yst / RstNorm;
float Rz = Zst / RstNorm;
float Xx = AntXaxisX / AntXaxisNorm;
float Xy = AntXaxisY / AntXaxisNorm;
float Xz = AntXaxisZ / AntXaxisNorm;
float Yx = AntYaxisX / AntYaxisNorm;
float Yy = AntYaxisY / AntYaxisNorm;
float Yz = AntYaxisZ / AntYaxisNorm;
float Zx = AntZaxisX / AntZaxisNorm;
float Zy = AntZaxisY / AntZaxisNorm;
float Zz = AntZaxisZ / AntZaxisNorm;
double Rx = Xst / RstNorm;
double Ry = Yst / RstNorm;
double Rz = Zst / RstNorm;
double Xx = AntXaxisX / AntXaxisNorm;
double Xy = AntXaxisY / AntXaxisNorm;
double Xz = AntXaxisZ / AntXaxisNorm;
double Yx = AntYaxisX / AntYaxisNorm;
double Yy = AntYaxisY / AntYaxisNorm;
double Yz = AntYaxisZ / AntYaxisNorm;
double Zx = AntZaxisX / AntZaxisNorm;
double Zy = AntZaxisY / AntZaxisNorm;
double Zz = AntZaxisZ / AntZaxisNorm;
float Xant = (Rx * Yy * Zz - Rx * Yz * Zy - Ry * Yx * Zz + Ry * Yz * Zx + Rz * Yx * Zy - Rz * Yy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
float Yant = -(Rx * Xy * Zz - Rx * Xz * Zy - Ry * Xx * Zz + Ry * Xz * Zx + Rz * Xx * Zy - Rz * Xy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
float Zant = (Rx * Xy * Yz - Rx * Xz * Yy - Ry * Xx * Yz + Ry * Xz * Yx + Rz * Xx * Yy - Rz * Xy * Yx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
double Xant = (Rx * Yy * Zz - Rx * Yz * Zy - Ry * Yx * Zz + Ry * Yz * Zx + Rz * Yx * Zy - Rz * Yy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
double Yant = -(Rx * Xy * Zz - Rx * Xz * Zy - Ry * Xx * Zz + Ry * Xz * Zx + Rz * Xx * Zy - Rz * Xy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
double Zant = (Rx * Xy * Yz - Rx * Xz * Yy - Ry * Xx * Yz + Ry * Xz * Yx + Rz * Xx * Yy - Rz * Xy * Yx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
// 计算theta 与 phi
float Norm = sqrtf(Xant * Xant + Yant * Yant + Zant * Zant); // 计算 pho
float ThetaAnt = acosf(Zant / Norm); // theta 与 Z轴的夹角
float PhiAnt = atanf(Yant / Xant); // -pi/2 ~pi/2
double Norm = sqrtf(Xant * Xant + Yant * Yant + Zant * Zant); // 计算 pho
double ThetaAnt = acosf(Zant / Norm); // theta 与 Z轴的夹角
double PhiAnt = atanf(Yant / Xant); // -pi/2 ~pi/2
if (abs(Yant) < PRECISIONTOLERANCE) { // X轴上
@ -278,17 +278,17 @@ __global__ void CUDA_SatelliteAntDirectNormal(float* RstX, float* RstY, float* R
}
}
__global__ void CUDA_BillerInterpAntPattern(float* antpattern,
float starttheta, float startphi, float dtheta, float dphi,
__global__ void CUDA_BillerInterpAntPattern(double* antpattern,
double starttheta, double startphi, double dtheta, double dphi,
long thetapoints, long phipoints,
float* searththeta, float* searchphi, float* searchantpattern,
double* searththeta, double* searchphi, double* searchantpattern,
long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float stheta = searththeta[idx];
float sphi = searchphi[idx];
float pthetaid = (stheta - starttheta) / dtheta;//
float pphiid = (sphi - startphi) / dphi;
double stheta = searththeta[idx];
double sphi = searchphi[idx];
double pthetaid = (stheta - starttheta) / dtheta;//
double pphiid = (sphi - startphi) / dphi;
long lasttheta = floorf(pthetaid);
long nextTheta = lasttheta + 1;
@ -301,18 +301,18 @@ __global__ void CUDA_BillerInterpAntPattern(float* antpattern,
searchantpattern[idx] = 0;
}
else {
float x = stheta;
float y = sphi;
double x = stheta;
double y = sphi;
float x1 = lasttheta * dtheta + starttheta;
float x2 = nextTheta * dtheta + starttheta;
float y1 = lastphi * dphi + startphi;
float y2 = nextPhi * dphi + startphi;
double x1 = lasttheta * dtheta + starttheta;
double x2 = nextTheta * dtheta + starttheta;
double y1 = lastphi * dphi + startphi;
double y2 = nextPhi * dphi + startphi;
float z11 = antpattern[lasttheta * phipoints + lastphi];
float z12 = antpattern[lasttheta * phipoints + nextPhi];
float z21 = antpattern[nextTheta * phipoints + lastphi];
float z22 = antpattern[nextTheta * phipoints + nextPhi];
double z11 = antpattern[lasttheta * phipoints + lastphi];
double z12 = antpattern[lasttheta * phipoints + nextPhi];
double z21 = antpattern[nextTheta * phipoints + lastphi];
double z22 = antpattern[nextTheta * phipoints + nextPhi];
z11 = powf(10, z11 / 10);
@ -320,7 +320,7 @@ __global__ void CUDA_BillerInterpAntPattern(float* antpattern,
z21 = powf(10, z21 / 10);
z22 = powf(10, z22 / 10);
float GainValue = (z11 * (x2 - x) * (y2 - y)
double GainValue = (z11 * (x2 - x) * (y2 - y)
+ z21 * (x - x1) * (y2 - y)
+ z12 * (x2 - x) * (y - y1)
+ z22 * (x - x1) * (y - y1));
@ -330,24 +330,24 @@ __global__ void CUDA_BillerInterpAntPattern(float* antpattern,
}
}
__global__ void CUDA_calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt,
float* localangle, float* R, float* slopeangle,
float nearRange, float Fs, float Pt, float lamda, long FreqIDmax,
__global__ void CUDA_calculationEcho(double* sigma0, double* TransAnt, double* ReciveAnt,
double* localangle, double* R, double* slopeangle,
double nearRange, double Fs, double Pt, double lamda, long FreqIDmax,
cuComplex* echoArr, long* FreqID,
long len) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float r = R[idx];
float amp = Pt * TransAnt[idx] * ReciveAnt[idx];
double r = R[idx];
double amp = Pt * TransAnt[idx] * ReciveAnt[idx];
amp = amp * sigma0[idx];
amp = amp / (powf(4 * LAMP_CUDA_PI, 2) * powf(r, 4)); // 反射强度
// 处理相位
float phi = (-4 * LAMP_CUDA_PI / lamda) * r;
double phi = (-4 * LAMP_CUDA_PI / lamda) * r;
cuComplex echophi = make_cuComplex(0, phi);
cuComplex echophiexp = cuCexpf(echophi);
float timeR = 2 * (r - nearRange) / LIGHTSPEED * Fs;
double timeR = 2 * (r - nearRange) / LIGHTSPEED * Fs;
long timeID = floorf(timeR);
//if (timeID < 0 || timeID >= FreqIDmax) {
// timeID = 0;
@ -361,15 +361,15 @@ __global__ void CUDA_calculationEcho(float* sigma0, float* TransAnt, float* Reci
}
__global__ void CUDA_AntPatternInterpGain(float* anttheta, float* antphi, float* gain,
float* antpattern, float starttheta, float startphi, float dtheta, float dphi, int thetapoints, int phipoints, long len) {
__global__ void CUDA_AntPatternInterpGain(double* anttheta, double* antphi, double* gain,
double* antpattern, double starttheta, double startphi, double dtheta, double dphi, int thetapoints, int phipoints, long len) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
float temptheta = anttheta[idx];
float tempphi = antphi[idx];
float antPatternGain = GPU_BillerInterpAntPattern(antpattern,
double temptheta = anttheta[idx];
double tempphi = antphi[idx];
double antPatternGain = GPU_BillerInterpAntPattern(antpattern,
starttheta, startphi, dtheta, dphi, thetapoints, phipoints,
temptheta, tempphi);
gain[idx] = antPatternGain;
@ -378,12 +378,12 @@ __global__ void CUDA_AntPatternInterpGain(float* anttheta, float* antphi, float*
__global__ void CUDA_InterpSigma(
long* demcls, float* sigmaAmp, float* localanglearr, long len,
long* demcls, double* sigmaAmp, double* localanglearr, long len,
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < len) {
long clsid = demcls[idx];
float localangle = localanglearr[idx];
double localangle = localanglearr[idx];
CUDASigmaParam tempsigma = sigma0Paramslist[clsid];
if (localangle < 0 || localangle >= LAMP_CUDA_PI / 2) {
sigmaAmp[idx] = 0;
@ -400,7 +400,7 @@ __global__ void CUDA_InterpSigma(
sigmaAmp[idx] = 0;
}
else {
float sigma = GPU_getSigma0dB(tempsigma, localangle);
double sigma = GPU_getSigma0dB(tempsigma, localangle);
sigma = powf(10.0, sigma / 10.0);// 后向散射系数
//printf("cls:%d;localangle=%f;sigma0=%f;\n", clsid, localangle, sigma);
sigmaAmp[idx] = sigma;
@ -413,23 +413,23 @@ __global__ void CUDA_InterpSigma(
__global__ void CUDAKernel_RFPC_Caluation_R_Gain(
float antX, float antY, float antZ, // 天线的坐标
float* targetX, float* targetY, float* targetZ, long len, // 地面坐标
double antX, double antY, double antZ, // 天线的坐标
double* targetX, double* targetY, double* targetZ, long len, // 地面坐标
long* demCls,
float* demSlopeX, float* demSlopeY, float* demSlopeZ, // 地表坡度矢量
float antXaxisX, float antXaxisY, float antXaxisZ, // 天线坐标系的X轴
float antYaxisX, float antYaxisY, float antYaxisZ,// 天线坐标系的Y轴
float antZaxisX, float antZaxisY, float antZaxisZ,// 天线坐标系的Z轴
float antDirectX, float antDirectY, float antDirectZ,// 天线的指向
float Pt,// 发射能量
double* demSlopeX, double* demSlopeY, double* demSlopeZ, // 地表坡度矢量
double antXaxisX, double antXaxisY, double antXaxisZ, // 天线坐标系的X轴
double antYaxisX, double antYaxisY, double antYaxisZ,// 天线坐标系的Y轴
double antZaxisX, double antZaxisY, double antZaxisZ,// 天线坐标系的Z轴
double antDirectX, double antDirectY, double antDirectZ,// 天线的指向
double Pt,// 发射能量
double refPhaseRange,
float* TransAntpattern, float Transtarttheta, float Transstartphi, float Transdtheta, float Transdphi, int Transthetapoints, int Transphipoints, // 发射天线方向图
float* ReceiveAntpattern, float Receivestarttheta, float Receivestartphi, float Receivedtheta, float Receivedphi, int Receivethetapoints, int Receivephipoints,//接收天线方向图
float NearR, float FarR, // 距离范围
double* TransAntpattern, double Transtarttheta, double Transstartphi, double Transdtheta, double Transdphi, int Transthetapoints, int Transphipoints, // 发射天线方向图
double* ReceiveAntpattern, double Receivestarttheta, double Receivestartphi, double Receivedtheta, double Receivedphi, int Receivethetapoints, int Receivephipoints,//接收天线方向图
double NearR, double FarR, // 距离范围
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 插值图
float* factorj, long freqnum,
double* factorj, long freqnum,
double* outR, // 输出距离
//float* outAmp // 输出增益
//double* outAmp // 输出增益
double* PRFEcho_real, double* PRFEcho_imag, long prfid
) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
@ -441,9 +441,9 @@ __global__ void CUDAKernel_RFPC_Caluation_R_Gain(
double RstY = antY - ty;
double RstZ = antZ - tz;
float slopeX = demSlopeX[idx];
float slopeY = demSlopeY[idx];
float slopeZ = demSlopeZ[idx];
double slopeX = demSlopeX[idx];
double slopeY = demSlopeY[idx];
double slopeZ = demSlopeZ[idx];
double RstR2 = RstX * RstX + RstY * RstY + RstZ * RstZ;
double RstR = sqrt(RstR2); // 矢量距离
@ -454,10 +454,10 @@ __global__ void CUDAKernel_RFPC_Caluation_R_Gain(
}
else {
// 求解坡度
float slopR = sqrtf(slopeX * slopeX + slopeY * slopeY + slopeZ * slopeZ); //
float dotAB = RstX * slopeX + RstY * slopeY + RstZ * slopeZ;
float localangle = acosf(dotAB / (RstR * slopR)); // 局地入射角
float ampGain = 0;
double slopR = sqrtf(slopeX * slopeX + slopeY * slopeY + slopeZ * slopeZ); //
double dotAB = RstX * slopeX + RstY * slopeY + RstZ * slopeZ;
double localangle = acosf(dotAB / (RstR * slopR)); // 局地入射角
double ampGain = 0;
// 求解天线方向图指向
CUDAVectorEllipsoidal antVector = GPU_SatelliteAntDirectNormal(
RstX, RstY, RstZ,
@ -468,22 +468,22 @@ __global__ void CUDAKernel_RFPC_Caluation_R_Gain(
);
if (antVector.Rho > 0) {
// 发射方向图
float temptheta = antVector.theta * r2d;
float tempphi = antVector.phi * r2d;
float TansantPatternGain =
double temptheta = antVector.theta * r2d;
double tempphi = antVector.phi * r2d;
double TansantPatternGain =
GPU_BillerInterpAntPattern(
TransAntpattern,
Transtarttheta, Transstartphi, Transdtheta, Transdphi, Transthetapoints, Transphipoints,
temptheta, tempphi);
// 接收方向图
float antPatternGain = GPU_BillerInterpAntPattern(
double antPatternGain = GPU_BillerInterpAntPattern(
ReceiveAntpattern,
Receivestarttheta, Receivestartphi, Receivedtheta, Receivedphi, Receivethetapoints, Receivephipoints,
temptheta, tempphi);
// 计算
float sigma0 = 0;
double sigma0 = 0;
{
long clsid = demCls[idx];
//printf("clsid=%d\n", clsid);
@ -503,25 +503,25 @@ __global__ void CUDAKernel_RFPC_Caluation_R_Gain(
sigma0 = 0;
}
else {
float sigma = GPU_getSigma0dB(tempsigma, localangle);
double sigma = GPU_getSigma0dB(tempsigma, localangle);
sigma0 = powf(10.0, sigma / 10.0);// 后向散射系数
}
}
ampGain = TansantPatternGain * antPatternGain;
ampGain = ampGain / (powf(4 * LAMP_CUDA_PI, 2) * powf(RstR, 4)); // 反射强度
float outAmp = ampGain * Pt * sigma0;
double outR = RstR- refPhaseRange;
double outAmp_temp = ampGain * Pt * sigma0;
double tempR = RstR- refPhaseRange;
outR[idx] = RstR ;
for (long ii = 0; ii < freqnum; ii++) {
float phi= outR * factorj[ii]; // 相位
double phi= tempR * factorj[ii]; // 相位
// Eular; exp(ix)=cos(x)+isin(x)
double real = outAmp * cos(phi); // 实部
double imag = outAmp * sin(phi); // 虚部
double real = outAmp_temp * cos(phi); // 实部
double imag = outAmp_temp * sin(phi); // 虚部
atomicAdd(&PRFEcho_real[prfid * freqnum+ ii], real);
atomicAdd(&PRFEcho_imag[prfid * freqnum+ ii], imag);
}
}
else {
}
@ -530,19 +530,19 @@ __global__ void CUDAKernel_RFPC_Caluation_R_Gain(
}
__global__ void CUDAKernel_PRF_CalFreqEcho(
double* Rarr, float* ampArr, long pixelcount,
float* factorj, long freqnum,
double* Rarr, double* ampArr, long pixelcount,
double* factorj, long freqnum,
double dx, double nearR,
cuComplex* PRFEcho, long prfid) {
long idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < freqnum) {
float fatorj = factorj[idx];
float phi = 0;
float amptemp = 0;
double fatorj = factorj[idx];
double phi = 0;
double amptemp = 0;
cuComplex tempfreqEcho = PRFEcho[prfid * freqnum + idx];
for (long i = 0; i < pixelcount; i++) { // 区域积分
//phi = (R = R - (floor(R / lamda) - 1) * lamda)* fatorj; // 相位
float phi = Rarr[i] * factorj[idx]; // 相位
double phi = Rarr[i] * factorj[idx]; // 相位
amptemp = ampArr[i];
//printf("amp=%f\n", amptemp);
// Eular; exp(ix)=cos(x)+isin(x)
@ -554,14 +554,73 @@ __global__ void CUDAKernel_PRF_CalFreqEcho(
}
}
__global__ void CUDAKernel_PRF_GeneratorEcho(double* Rarr, double* ampArr, long blocknum, long pixelcount, double* factorj, long freqnum,
double nearR, double farR, double* echo_real, double* echo_imag, long prfid) //11
{
//// 假定共享内存大小为49152 byte
//// 假定每个Block 线程数大小为 256
__shared__ double s_R[GPU_SHARE_MEMORY]; // 距离 256*12 * 8= 49.2kb
__shared__ double s_Amp[GPU_SHARE_MEMORY]; // 振幅 3072 * 8= 49.2kb 49.2*2 = 98.4 < 100 KB
const int bid = blockIdx.x; // 获取 grid网格编号ID
const int tid = threadIdx.x;// 获取 单个 block 中的线程ID
const int startPIX = bid * GPU_SHARE_MEMORY;
int curthreadidx = 0;
for (long i = 0; i < GPU_SHARE_STEP; i++) {
curthreadidx = tid * GPU_SHARE_STEP + i;
s_R[curthreadidx] = (startPIX + curthreadidx) < pixelcount ? Rarr[startPIX + curthreadidx] : 0.0;
s_Amp[curthreadidx] = (startPIX + curthreadidx) < pixelcount ? ampArr[startPIX + curthreadidx] : 0.0;
}
__syncthreads(); // 确定所有待处理数据都已经进入程序中
long freqnumblock = freqnum / 256 + 1; //16
//if (startPIX < pixelcount) { // 存在可能处理的计算
// double temp_real = 0;
// double temp_imag = 0;
// double factorjTemp = 0;
// double temp_phi = 0;
// double temp_amp = 0;
// curthreadidx = 0;
// for (long i = 0; i < freqnumblock; i++) {
// curthreadidx = tid * freqnumblock + i; // 获取当前频率
// if (curthreadidx < freqnum) { // 存在频率
// factorjTemp = factorj[curthreadidx];
// for (long j = 0; j < GPU_SHARE_MEMORY; j++) {
// temp_phi = s_R[j] * factorjTemp;
// temp_amp = s_Amp[j];
// temp_real = temp_real + temp_amp * cos(temp_phi);
// temp_imag = temp_imag + temp_amp * sin(temp_phi);
// }
// //atomicAdd(&echo_real[prfid * freqnum + curthreadidx], temp_real); // 更新实部
// //atomicAdd(&echo_imag[prfid * freqnum + curthreadidx], temp_imag); // 更新虚部
// }
// }
//}
}
/** 对外封装接口 *******************************************************************************************************/
extern "C" void SatelliteAntDirectNormal(float* RstX, float* RstY, float* RstZ,
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ,
float* thetaAnt, float* phiAnt,
extern "C" void SatelliteAntDirectNormal(double* RstX, double* RstY, double* RstZ,
double antXaxisX, double antXaxisY, double antXaxisZ,
double antYaxisX, double antYaxisY, double antYaxisZ,
double antZaxisX, double antZaxisY, double antZaxisZ,
double antDirectX, double antDirectY, double antDirectZ,
double* thetaAnt, double* phiAnt,
long len) {
int blockSize = 256; // 每个块的线程数
@ -586,8 +645,8 @@ extern "C" void SatelliteAntDirectNormal(float* RstX, float* RstY, float* RstZ,
}
extern "C" void AntPatternInterpGain(float* anttheta, float* antphi, float* gain,
float* antpattern, float starttheta, float startphi, float dtheta, float dphi, int thetapoints, int phipoints, long len) {
extern "C" void AntPatternInterpGain(double* anttheta, double* antphi, double* gain,
double* antpattern, double starttheta, double startphi, double dtheta, double dphi, int thetapoints, int phipoints, long len) {
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
//printf("\nCUDA_RFPC_SiglePRF blockSize:%d ,numBlock:%d\n", blockSize, numBlocks);
@ -607,9 +666,9 @@ extern "C" void AntPatternInterpGain(float* anttheta, float* antphi, float* ga
}
extern "C" void calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt,
float* localangle, float* R, float* slopeangle,
float nearRange, float Fs, float pt, float lamda, long FreqIDmax,
extern "C" void calculationEcho(double* sigma0, double* TransAnt, double* ReciveAnt,
double* localangle, double* R, double* slopeangle,
double nearRange, double Fs, double pt, double lamda, long FreqIDmax,
cuComplex* echoAmp, long* FreqID,
long len)
{
@ -635,7 +694,7 @@ extern "C" void calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt
extern "C" void CUDAInterpSigma(
long* demcls, float* sigmaAmp, float* localanglearr, long len,
long* demcls, double* sigmaAmp, double* localanglearr, long len,
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen) {// 地表覆盖类型-sigma插值对应函数-ulaby
int blockSize = 256; // 每个块的线程数
int numBlocks = (len + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
@ -655,23 +714,23 @@ extern "C" void CUDAInterpSigma(
}
extern "C" void CUDARFPC_Caluation_R_Gain(
float antX, float antY, float antZ, // 天线的坐标
float* targetX, float* targetY, float* targetZ, long TargetPixelNumber, // 地面坐标
double antX, double antY, double antZ, // 天线的坐标
double* targetX, double* targetY, double* targetZ, long TargetPixelNumber, // 地面坐标
long* demCls,
float* demSlopeX, float* demSlopeY, float* demSlopeZ, // 地表坡度矢量
float antXaxisX, float antXaxisY, float antXaxisZ, // 天线坐标系的X轴
float antYaxisX, float antYaxisY, float antYaxisZ,// 天线坐标系的Y轴
float antZaxisX, float antZaxisY, float antZaxisZ,// 天线坐标系的Z轴
float antDirectX, float antDirectY, float antDirectZ,// 天线的指向
float Pt,// 发射能量
double* demSlopeX, double* demSlopeY, double* demSlopeZ, // 地表坡度矢量
double antXaxisX, double antXaxisY, double antXaxisZ, // 天线坐标系的X轴
double antYaxisX, double antYaxisY, double antYaxisZ,// 天线坐标系的Y轴
double antZaxisX, double antZaxisY, double antZaxisZ,// 天线坐标系的Z轴
double antDirectX, double antDirectY, double antDirectZ,// 天线的指向
double Pt,// 发射能量
double refPhaseRange,
float* TransAntpattern, float Transtarttheta, float Transstartphi, float Transdtheta, float Transdphi, int Transthetapoints, int Transphipoints, // 发射天线方向图
float* ReceiveAntpattern, float Receivestarttheta, float Receivestartphi, float Receivedtheta, float Receivedphi, int Receivethetapoints, int Receivephipoints,//接收天线方向图
float NearR, float FarR, // 距离范围
double* TransAntpattern, double Transtarttheta, double Transstartphi, double Transdtheta, double Transdphi, int Transthetapoints, int Transphipoints, // 发射天线方向图
double* ReceiveAntpattern, double Receivestarttheta, double Receivestartphi, double Receivedtheta, double Receivedphi, int Receivethetapoints, int Receivephipoints,//接收天线方向图
double NearR, double FarR, // 距离范围
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 插值图
float* factorj, long freqnum,
double* factorj, long freqnum,
double* outR, // 输出距离
//float* outAmp // 输出增益
//double* outAmp // 输出增益
double* PRFEcho_real, double* PRFEcho_imag, long prfid
)
{
@ -714,8 +773,8 @@ extern "C" void CUDARFPC_Caluation_R_Gain(
}
extern "C" void CUDA_PRF_CalFreqEcho(
double* Rarr, float* ampArr, long pixelcount,
float* factorj, long freqnum,
double* Rarr, double* ampArr, long pixelcount,
double* factorj, long freqnum,
double dx, double nearR,
cuComplex* PRFEcho, long prfid)
{
@ -738,6 +797,46 @@ extern "C" void CUDA_PRF_CalFreqEcho(
cudaDeviceSynchronize();
}
extern "C" void CUDA_PRF_GeneratorEcho(cublasHandle_t handle,double* Rarr, double* ampArr, long pixelcount, double* factorj, long freqnum, double nearR, double farR, double* echo_real, double* echo_imag, long prfid)
{
//cublasHandle_t handle;
//cublasStatus_t status = cublasCreate(&handle);
long blocknum = pixelcount / GPU_SHARE_MEMORY + 1;
int blockSize = 256; // 每个块的线程数
int numBlocks = (pixelcount + GPU_SHARE_MEMORY - 1) / GPU_SHARE_MEMORY; // 网格数量
CUDAKernel_PRF_GeneratorEcho << <numBlocks, blockSize >> > (Rarr, ampArr, blocknum, pixelcount,
factorj, freqnum,
nearR, farR,
echo_real, echo_imag, prfid);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDA_PRF_GeneratorEcho CUDA Error: %s\n", cudaGetErrorString(err));
// Possibly: exit(-1) if program cannot continue....
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
//cublasDestroy(handle);
}
void CUDA_RFPC_MainProgramm()
{
// 创建 cuBLAS 句柄
cublasHandle_t handle;
cublasStatus_t status = cublasCreate(&handle);
cublasDestroy(handle);
}
#endif

View File

@ -10,162 +10,145 @@
extern "C" struct CUDASigmaParam {
float p1;
float p2;
float p3;
float p4;
float p5;
float p6;
double p1;
double p2;
double p3;
double p4;
double p5;
double p6;
};
extern __device__ float GPU_getSigma0dB(CUDASigmaParam param, float theta);
extern __device__ double GPU_getSigma0dB(CUDASigmaParam param, double theta);
extern __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
float RstX, float RstY, float RstZ,
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ
double RstX, double RstY, double RstZ,
double antXaxisX, double antXaxisY, double antXaxisZ,
double antYaxisX, double antYaxisY, double antYaxisZ,
double antZaxisX, double antZaxisY, double antZaxisZ,
double antDirectX, double antDirectY, double antDirectZ
);
extern __device__ float GPU_BillerInterpAntPattern(float* antpattern,
float starttheta, float startphi, float dtheta, float dphi,
extern __device__ double GPU_BillerInterpAntPattern(double* antpattern,
double starttheta, double startphi, double dtheta, double dphi,
long thetapoints, long phipoints,
float searththeta, float searchphi);
double searththeta, double searchphi);
extern __global__ void CUDA_AntPatternInterpGain(float* anttheta, float* antphi, float* gain,
float* antpattern, float starttheta, float startphi, float dtheta, float dphi, int thetapoints, int phipoints, long len);
extern __global__ void CUDA_AntPatternInterpGain(double* anttheta, double* antphi, double* gain,
double* antpattern, double starttheta, double startphi, double dtheta, double dphi, int thetapoints, int phipoints, long len);
extern __global__ void CUDA_InterpSigma(
long* demcls, float* sigmaAmp, float* localanglearr, long len,
long* demcls, double* sigmaAmp, double* localanglearr, long len,
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen);
extern __global__ void CUDAKernel_RFPC_Caluation_R_Gain(
float antX, float antY, float antZ, // 天线的坐标
float* targetX, float* targetY, float* targetZ, long len, // 地面坐标
long* demCls,
float* demSlopeX, float* demSlopeY, float* demSlopeZ, // 地表坡度矢量
float antXaxisX, float antXaxisY, float antXaxisZ, // 天线坐标系的X轴
float antYaxisX, float antYaxisY, float antYaxisZ,// 天线坐标系的Y轴
float antZaxisX, float antZaxisY, float antZaxisZ,// 天线坐标系的Z轴
float antDirectX, float antDirectY, float antDirectZ,// 天线的指向
float Pt,// 发射能量
double refPhaseRange,
float* TransAntpattern, float Transtarttheta, float Transstartphi, float Transdtheta, float Transdphi, int Transthetapoints, int Transphipoints, // 发射天线方向图
float* ReceiveAntpattern, float Receivestarttheta, float Receivestartphi, float Receivedtheta, float Receivedphi, int Receivethetapoints, int Receivephipoints,//接收天线方向图
float NearR, float FarR, // 距离范围
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 插值图
float* factorj, long freqnum,
double* outR, // 输出距离
float* outAmp // 输出增益
);
extern __global__ void CUDAKernel_PRF_CalFreqEcho(
double* Rarr, float* ampArr, long pixelcount,
float* factorj, long freqnum,
double* Rarr, double* ampArr, long pixelcount,
double* factorj, long freqnum,
double dx, double nearR,
cuComplex* PRFEcho, long prfid);
extern __device__ cuComplex GPU_calculationEcho(float sigma0, float TransAnt, float ReciveAnt,
float localangle, float R, float slopeangle, float Pt, float lamda);
extern __device__ cuComplex GPU_calculationEcho(double sigma0, double TransAnt, double ReciveAnt,
double localangle, double R, double slopeangle, double Pt, double lamda);
extern __global__ void CUDA_calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt,
float* localangle, float* R, float* slopeangle,
float nearRange, float Fs, float Pt, float lamda, long FreqIDmax,
extern __global__ void CUDA_calculationEcho(double* sigma0, double* TransAnt, double* ReciveAnt,
double* localangle, double* R, double* slopeangle,
double nearRange, double Fs, double Pt, double lamda, long FreqIDmax,
cuComplex* echoArr, long* FreqID,
long len);
extern __global__ void CUDA_BillerInterpAntPattern(float* antpattern,
float starttheta, float startphi, float dtheta, float dphi,
extern __global__ void CUDA_BillerInterpAntPattern(double* antpattern,
double starttheta, double startphi, double dtheta, double dphi,
long thetapoints, long phipoints,
float* searththeta, float* searchphi, float* searchantpattern,
double* searththeta, double* searchphi, double* searchantpattern,
long len);
extern __global__ void CUDA_SatelliteAntDirectNormal(float* RstX, float* RstY, float* RstZ,
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ,
float* thetaAnt, float* phiAnt
extern __global__ void CUDA_SatelliteAntDirectNormal(double* RstX, double* RstY, double* RstZ,
double antXaxisX, double antXaxisY, double antXaxisZ,
double antYaxisX, double antYaxisY, double antYaxisZ,
double antZaxisX, double antZaxisY, double antZaxisZ,
double antDirectX, double antDirectY, double antDirectZ,
double* thetaAnt, double* phiAnt
, long len);
extern __global__ void CUDAKernel_PRFSumEcho_Rows(
double* Rarr, float* ampArr, long Rows, long Cols,
long startRid,
float* factorj, long freqnum,
cuComplex* freqRowsbuffer, long tempRows
);
extern __global__ void CUDAKernel_PRF_GeneratorEcho(double* Rarr, double* ampArr, long blocknum, long pixelcount, double* factorj, long freqnum,
double nearR, double farR, double* echo_real, double* echo_imag, long prfid);
extern __global__ void CUDAKernel_PRFSumEcho_Freq(
cuComplex* freqRowsbuffer, long tempRows, long freqnum,
cuComplex* PRFEcho, long prfid
);
extern "C" void SatelliteAntDirectNormal(double* RstX, double* RstY, double* RstZ,
double antXaxisX, double antXaxisY, double antXaxisZ,
double antYaxisX, double antYaxisY, double antYaxisZ,
double antZaxisX, double antZaxisY, double antZaxisZ,
double antDirectX, double antDirectY, double antDirectZ,
double* thetaAnt, double* phiAnt, long len);
extern "C" void SatelliteAntDirectNormal(float* RstX, float* RstY, float* RstZ,
float antXaxisX, float antXaxisY, float antXaxisZ,
float antYaxisX, float antYaxisY, float antYaxisZ,
float antZaxisX, float antZaxisY, float antZaxisZ,
float antDirectX, float antDirectY, float antDirectZ,
float* thetaAnt, float* phiAnt, long len);
extern "C" void AntPatternInterpGain(float* anttheta, float* antphi, float* gain,
float* antpattern,
float starttheta, float startphi, float dtheta, float dphi, int thetapoints, int phipoints,
extern "C" void AntPatternInterpGain(double* anttheta, double* antphi, double* gain,
double* antpattern,
double starttheta, double startphi, double dtheta, double dphi, int thetapoints, int phipoints,
long len);
extern "C" void calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt,
float* localangle, float* R, float* slopeangle,
float nearRange, float Fs, float pt, float lamda, long FreqIDmax,
extern "C" void calculationEcho(double* sigma0, double* TransAnt, double* ReciveAnt,
double* localangle, double* R, double* slopeangle,
double nearRange, double Fs, double pt, double lamda, long FreqIDmax,
cuComplex* echoAmp, long* FreqID,
long len);
extern "C" void CUDAInterpSigma(
long* demcls, float* sigmaAmp, float* localanglearr, long len,
long* demcls, double* sigmaAmp, double* localanglearr, long len,
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen);
// 计算坐标的 距离、增益
extern "C" void CUDARFPC_Caluation_R_Gain(
float antX, float antY, float antZ, // 天线的坐标
float* targetX, float* targetY, float* targetZ, long TargetPixelNumber, // 地面坐标
double antX, double antY, double antZ, // 天线的坐标
double* targetX, double* targetY, double* targetZ, long TargetPixelNumber, // 地面坐标
long* demCls,
float* demSlopeX, float* demSlopeY, float* demSlopeZ, // 地表坡度矢量
float antXaxisX, float antXaxisY, float antXaxisZ, // 天线坐标系的X轴
float antYaxisX, float antYaxisY, float antYaxisZ,// 天线坐标系的Y轴
float antZaxisX, float antZaxisY, float antZaxisZ,// 天线坐标系的Z轴
float antDirectX, float antDirectY, float antDirectZ,// 天线的指向
float Pt,// 发射能量
double* demSlopeX, double* demSlopeY, double* demSlopeZ, // 地表坡度矢量
double antXaxisX, double antXaxisY, double antXaxisZ, // 天线坐标系的X轴
double antYaxisX, double antYaxisY, double antYaxisZ,// 天线坐标系的Y轴
double antZaxisX, double antZaxisY, double antZaxisZ,// 天线坐标系的Z轴
double antDirectX, double antDirectY, double antDirectZ,// 天线的指向
double Pt,// 发射能量
double refPhaseRange,
float* TransAntpattern, float Transtarttheta, float Transstartphi, float Transdtheta, float Transdphi, int Transthetapoints, int Transphipoints, // 发射天线方向图
float* ReceiveAntpattern, float Receivestarttheta, float Receivestartphi, float Receivedtheta, float Receivedphi, int Receivethetapoints, int Receivephipoints,//接收天线方向图
float NearR, float FarR, // 距离范围
double* TransAntpattern, double Transtarttheta, double Transstartphi, double Transdtheta, double Transdphi, int Transthetapoints, int Transphipoints, // 发射天线方向图
double* ReceiveAntpattern, double Receivestarttheta, double Receivestartphi, double Receivedtheta, double Receivedphi, int Receivethetapoints, int Receivephipoints,//接收天线方向图
double NearR, double FarR, // 距离范围
CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 插值图
float* factorj, long freqnum,
double* factorj, long freqnum,
double* outR, // 输出距离
//float* outAmp // 输出增益
//double* outAmp // 输出增益
double* PRFEcho_real, double* PRFEcho_imag, long prfid
);
extern "C" void CUDA_PRF_CalFreqEcho(
double* Rarr, float* ampArr, long pixelcount,
float* factorj, long freqnum,
double* Rarr, double* ampArr, long pixelcount,
double* factorj, long freqnum,
double dx, double nearR,
cuComplex* PRFEcho, long prfid);
extern "C" void CUDA_PRF_GeneratorEcho(cublasHandle_t handle,double* Rarr, double* ampArr, long pixelcount,
double* factorj, long freqnum,
double nearR, double farR,
double* echo_real, double* echo_imag, long prfid);
#endif

View File

@ -15,6 +15,8 @@
#define BLOCK_SIZE 256
#define GPU_SHARE_MEMORY 2816
#define GPU_SHARE_STEP 11
// 默认显存分布

View File

@ -69,6 +69,7 @@
</PropertyGroup>
<PropertyGroup Condition="'$(Configuration)|$(Platform)' == 'Release|x64'">
<IncludePath>.\LAMPScatterTool;.\GPUTool;.\SimulationSAR;.\GF3ProcessToolbox;.\BaseTool;$(oneMKLIncludeDir);$(IncludePath)</IncludePath>
<ReferencePath>C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.6\lib\x64;$(ReferencePath)</ReferencePath>
</PropertyGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'">
<ClCompile>
@ -83,9 +84,11 @@
<Link>
<LinkTimeCodeGeneration>Default</LinkTimeCodeGeneration>
<LargeAddressAware>true</LargeAddressAware>
<AdditionalDependencies>cublas.lib;%(AdditionalDependencies)</AdditionalDependencies>
</Link>
<CudaCompile>
<CodeGeneration>compute_86,sm_86</CodeGeneration>
<GenerateRelocatableDeviceCode>true</GenerateRelocatableDeviceCode>
</CudaCompile>
</ItemDefinitionGroup>
<ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)' == 'Debug|x64'" Label="Configuration">

View File

@ -22,6 +22,8 @@
#ifdef __CUDANVCC___
#include "GPUTool.cuh"
#include "GPURFPC.cuh"
#include <cuda_runtime.h>
#include <cublas_v2.h>
#endif // __CUDANVCC___
#include <Imageshow/ImageShowDialogClass.h>
@ -200,7 +202,7 @@ ErrorCode RFPCProcessCls::DEMPreprocess()
this->demxyzPath = QDir(tmpfolderPath).filePath("demxyz.bin");
gdalImage demds(this->DEMTiffPath);
gdalImage demxyz = CreategdalImage(demxyzPath, demds.height, demds.width, 3, demds.gt, demds.projection, true, true, true);// X,Y,Z
gdalImage demxyz = CreategdalImageDouble(demxyzPath, demds.height, demds.width, 3, demds.gt, demds.projection, true, true, true);// X,Y,Z
// 分块计算并转换为XYZ
@ -247,7 +249,7 @@ ErrorCode RFPCProcessCls::DEMPreprocess()
this->demsloperPath = QDir(tmpfolderPath).filePath("demsloper.bin");
this->demmaskPath = QDir(tmpfolderPath).filePath("demmask.bin");
gdalImage demsloperxyz = CreategdalImage(this->demsloperPath, demds.height, demds.width, 4, demds.gt, demds.projection, true, true, true);// X,Y,Z,cosangle
gdalImage demsloperxyz = CreategdalImageDouble(this->demsloperPath, demds.height, demds.width, 4, demds.gt, demds.projection, true, true, true);// X,Y,Z,cosangle
gdalImage demmask = CreategdalImage(this->demmaskPath, demxyz.height, demxyz.width, 1, demxyz.gt, demxyz.projection, true, true, true);// X,Y,Z
line_invert = 1000;
@ -350,6 +352,14 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
double R = 0;
double dem_row = 0, dem_col = 0, dem_alt = 0;
QVector<double> freqlist = this->TaskSetting->getFreqList();
long freqnum = freqlist.count();
std::shared_ptr<double> freqPtr(new double[freqnum], delArrPtr);
for (long ii = 0; ii < freqlist.count(); ii++) {
freqPtr.get()[ii] = freqlist[ii];
}
testOutAmpArr("freqlist.bin", (double*)(freqPtr.get()), freqnum, 1);
long double imageStarttime = 0;
imageStarttime = this->TaskSetting->getSARImageStartTime();
@ -359,17 +369,17 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
// 回波
long echoIdx = 0;
float NearRange = this->EchoSimulationData->getNearRange(); // 近斜距
float FarRange = this->EchoSimulationData->getFarRange();
double NearRange = this->EchoSimulationData->getNearRange(); // 近斜距
double FarRange = this->EchoSimulationData->getFarRange();
float TimgNearRange = 2 * NearRange / LIGHTSPEED;
float TimgFarRange = 2 * FarRange / LIGHTSPEED;
float dx = (FarRange - NearRange) / (PlusePoint - 1);
float Fs = this->TaskSetting->getFs(); // 距离向采样率
float Pt = this->TaskSetting->getPt() * this->TaskSetting->getGri();// 发射电压 1v
double TimgNearRange = 2 * NearRange / LIGHTSPEED;
double TimgFarRange = 2 * FarRange / LIGHTSPEED;
double dx = (FarRange - NearRange) / (PlusePoint - 1);
double Fs = this->TaskSetting->getFs(); // 距离向采样率
double Pt = this->TaskSetting->getPt() * this->TaskSetting->getGri();// 发射电压 1v
//double GainAntLen = -3;// -3dB 为天线半径
long pluseCount = this->PluseCount;
float lamda = this->TaskSetting->getCenterLamda(); // 波长
double lamda = this->TaskSetting->getCenterLamda(); // 波长
double refphaseRange = this->TaskSetting->getRefphaseRange(); // 参考相位斜距
// 天线方向图
std::shared_ptr<AbstractRadiationPattern> TransformPattern = this->TaskSetting->getTransformRadiationPattern(); // 发射天线方向图
@ -419,8 +429,8 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
double Tdtheta = (Tmaxtheta - Tmintheta) / (Tthetanum - 1);
double Tdphi = (Tmaxphi - Tminphi) / (Tphinum - 1);
float* h_TantPattern = (float*)mallocCUDAHost(sizeof(float) * Tthetanum * Tphinum);
float* d_TantPattern = (float*)mallocCUDADevice(sizeof(float) * Tthetanum * Tphinum);
double* h_TantPattern = (double*)mallocCUDAHost(sizeof(double) * Tthetanum * Tphinum);
double* d_TantPattern = (double*)mallocCUDADevice(sizeof(double) * Tthetanum * Tphinum);
for (long i = 0; i < Tthetanum; i++) {
for (long j = Tphinum - 1; j >= 0; j--) {
@ -435,7 +445,7 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
h_TantPattern[i * Tphinum + j] = powf(10.0, h_TantPattern[i * Tphinum + j] / 10);
}
}
HostToDevice(h_TantPattern, d_TantPattern, sizeof(float) * Tthetanum * Tphinum);
HostToDevice(h_TantPattern, d_TantPattern, sizeof(double) * Tthetanum * Tphinum);
// 处理接收天线方向图
double Rminphi = ReceivePattern->getMinPhi();
@ -452,8 +462,8 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
double Rdtheta = (Rmaxtheta - Rmintheta) / (Rthetanum - 1);
double Rdphi = (Rmaxphi - Rminphi) / (Rphinum - 1);
float* h_RantPattern = (float*)mallocCUDAHost(sizeof(float) * Rthetanum * Rphinum);
float* d_RantPattern = (float*)mallocCUDADevice(sizeof(float) * Rthetanum * Rphinum);
double* h_RantPattern = (double*)mallocCUDAHost(sizeof(double) * Rthetanum * Rphinum);
double* d_RantPattern = (double*)mallocCUDADevice(sizeof(double) * Rthetanum * Rphinum);
for (long i = 0; i < Rthetanum; i++) {
for (long j = 0; j < Rphinum; j++) {
@ -469,7 +479,7 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
h_RantPattern[i * Tphinum + j] = powf(10.0, h_RantPattern[i * Tphinum + j] / 10);
}
}
HostToDevice(h_RantPattern, d_RantPattern, sizeof(float) * Rthetanum * Rphinum);
HostToDevice(h_RantPattern, d_RantPattern, sizeof(double) * Rthetanum * Rphinum);
//处理地表覆盖
QMap<long, long> clamap;
@ -547,32 +557,44 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
Eigen::MatrixXd demsloper_y = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
Eigen::MatrixXd demsloper_z = Eigen::MatrixXd::Zero(tempDemRows, tempDemCols);
float* h_dem_x = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_dem_y = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_dem_z = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_demsloper_x = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_demsloper_y = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* h_demsloper_z = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
double* h_dem_x = (double*)mallocCUDAHost(sizeof(double) * blokline * tempDemCols);
double* h_dem_y = (double*)mallocCUDAHost(sizeof(double) * blokline * tempDemCols);
double* h_dem_z = (double*)mallocCUDAHost(sizeof(double) * blokline * tempDemCols);
double* h_demsloper_x = (double*)mallocCUDAHost(sizeof(double) * blokline * tempDemCols);
double* h_demsloper_y = (double*)mallocCUDAHost(sizeof(double) * blokline * tempDemCols);
double* h_demsloper_z = (double*)mallocCUDAHost(sizeof(double) * blokline * tempDemCols);
float* d_dem_x = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols); // 7
float* d_dem_y = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
float* d_dem_z = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
float* d_demsloper_x = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
float* d_demsloper_y = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
float* d_demsloper_z = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
double* d_dem_x = (double*)mallocCUDADevice(sizeof(double) * blokline * tempDemCols); // 7
double* d_dem_y = (double*)mallocCUDADevice(sizeof(double) * blokline * tempDemCols);
double* d_dem_z = (double*)mallocCUDADevice(sizeof(double) * blokline * tempDemCols);
double* d_demsloper_x = (double*)mallocCUDADevice(sizeof(double) * blokline * tempDemCols);
double* d_demsloper_y = (double*)mallocCUDADevice(sizeof(double) * blokline * tempDemCols);
double* d_demsloper_z = (double*)mallocCUDADevice(sizeof(double) * blokline * tempDemCols);
// 提前声明参数变量
double* h_R = (double*)mallocCUDAHost(sizeof(double) * blokline * tempDemCols);
double* d_R = (double*)mallocCUDADevice(sizeof(double) * blokline * tempDemCols);
float* h_amp = (float*)mallocCUDAHost(sizeof(float) * blokline * tempDemCols);
float* d_amp = (float*)mallocCUDADevice(sizeof(float) * blokline * tempDemCols);
double* h_amp = (double*)mallocCUDAHost(sizeof(double) * blokline * tempDemCols);
double* d_amp = (double*)mallocCUDADevice(sizeof(double) * blokline * tempDemCols);
double* h_one = (double*)mallocCUDAHost(sizeof(double) * blokline * tempDemCols);
double* d_one = (double*)mallocCUDADevice(sizeof(double) * blokline * tempDemCols);
for (long ii = 0; ii < blokline * tempDemCols; ii++) {
h_one[ii] = 1;
}
HostToDevice(h_one, d_one, sizeof(double) * blokline * tempDemCols);
double* h_temp = (double*)mallocCUDAHost(sizeof(double) * blokline * tempDemCols);
double* d_temp = (double*)mallocCUDADevice(sizeof(double) * blokline * tempDemCols);
// 地面回波
cuComplex* h_echo = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * blokline * tempDemCols);
cuComplex* d_echo = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * blokline * tempDemCols); //19
long echoblockline = Memory1GB / 8 / 2 / PlusePoint * 2;
long echoblockline = Memory1GB / 8 / 2 / PlusePoint * 1;
// 每一行的脉冲
double* h_PRFEcho_real = (double*)mallocCUDAHost(sizeof(double) * echoblockline * PlusePoint);
@ -580,18 +602,23 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
double* d_PRFEcho_real = (double*)mallocCUDADevice(sizeof(double) * echoblockline * PlusePoint);
double* d_PRFEcho_imag = (double*)mallocCUDADevice(sizeof(double) * echoblockline * PlusePoint);
float* h_factorj = (float*)mallocCUDAHost(sizeof(float) * freqlist.size());
float* d_factorj = (float*)mallocCUDADevice(sizeof(float) * freqlist.size());
double* h_factorj = (double*)mallocCUDAHost(sizeof(double) * freqlist.size());
double* d_factorj = (double*)mallocCUDADevice(sizeof(double) * freqlist.size());
for (long ii = 0; ii < freqlist.size(); ii++) {
h_factorj[ii] = -4*PI*freqlist[ii]/LIGHTSPEED;
}
HostToDevice(h_factorj, d_factorj, sizeof(float) * freqlist.size());
HostToDevice(h_factorj, d_factorj, sizeof(double) * freqlist.size());
testOutAmpArr("factorj.bin", h_factorj, freqlist.size(), 1);
// 地表覆盖类型
Eigen::MatrixXd landcover = Eigen::MatrixXd::Zero(blokline, tempDemCols);// 地面覆盖类型
long* h_demcls = (long*)mallocCUDAHost(sizeof(long) * blokline * tempDemCols);
long* d_demcls = (long*)mallocCUDADevice(sizeof(long) * blokline * tempDemCols);
cublasHandle_t handle;
cublasStatus_t status = cublasCreate(&handle);
for (startline = 0; startline < demRow; startline = startline + blokline) {
long newblokline = blokline;
if ((startline + blokline) >= demRow) {
@ -627,58 +654,70 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
FreeCUDAHost(h_amp); FreeCUDADevice(d_amp);
FreeCUDAHost(h_echo); FreeCUDADevice(d_echo);//19
FreeCUDAHost(h_demcls); FreeCUDADevice(d_demcls);
FreeCUDAHost(h_one); FreeCUDADevice(d_one);
FreeCUDAHost(h_temp); FreeCUDADevice(d_temp);
h_dem_x = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_dem_y = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_dem_z = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_demsloper_x = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_demsloper_y = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_demsloper_z = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_dem_x = (double*)mallocCUDAHost(sizeof(double) * newblokline * tempDemCols);
h_dem_y = (double*)mallocCUDAHost(sizeof(double) * newblokline * tempDemCols);
h_dem_z = (double*)mallocCUDAHost(sizeof(double) * newblokline * tempDemCols);
h_demsloper_x = (double*)mallocCUDAHost(sizeof(double) * newblokline * tempDemCols);
h_demsloper_y = (double*)mallocCUDAHost(sizeof(double) * newblokline * tempDemCols);
h_demsloper_z = (double*)mallocCUDAHost(sizeof(double) * newblokline * tempDemCols);
h_R = (double*)mallocCUDAHost(sizeof(double) * newblokline * tempDemCols);
h_amp = (float*)mallocCUDAHost(sizeof(float) * newblokline * tempDemCols);
h_amp = (double*)mallocCUDAHost(sizeof(double) * newblokline * tempDemCols);
h_echo = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * newblokline * tempDemCols);
h_demcls = (long*)mallocCUDAHost(sizeof(long) * newblokline * tempDemCols);
d_dem_x = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_dem_y = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_dem_z = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demsloper_x = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demsloper_y = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_demsloper_z = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);//6
d_amp = (float*)mallocCUDADevice(sizeof(float) * newblokline * tempDemCols);
d_dem_x = (double*)mallocCUDADevice(sizeof(double) * newblokline * tempDemCols);
d_dem_y = (double*)mallocCUDADevice(sizeof(double) * newblokline * tempDemCols);
d_dem_z = (double*)mallocCUDADevice(sizeof(double) * newblokline * tempDemCols);
d_demsloper_x = (double*)mallocCUDADevice(sizeof(double) * newblokline * tempDemCols);
d_demsloper_y = (double*)mallocCUDADevice(sizeof(double) * newblokline * tempDemCols);
d_demsloper_z = (double*)mallocCUDADevice(sizeof(double) * newblokline * tempDemCols);//6
d_amp = (double*)mallocCUDADevice(sizeof(double) * newblokline * tempDemCols);
d_echo = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * newblokline * tempDemCols);
d_demcls = (long*)mallocCUDADevice(sizeof(long) * newblokline * tempDemCols);
d_R = (double*)mallocCUDADevice(sizeof(double) * newblokline * tempDemCols);
h_one = (double*)mallocCUDAHost(sizeof(double) * newblokline * tempDemCols);
d_one = (double*)mallocCUDADevice(sizeof(double) * newblokline * tempDemCols);
h_temp = (double*)mallocCUDAHost(sizeof(double) * newblokline * tempDemCols);
d_temp = (double*)mallocCUDADevice(sizeof(double) * newblokline * tempDemCols);
}
//# pragma omp parallel for
for (long i = 0; i < newblokline; i++) {
for (long j = 0; j < blockwidth; j++) {
#ifdef __PRFDEBUG__
h_dem_x[i * blockwidth + j] = -2028380.6250000; float(dem_x(i, j));
h_dem_y[i * blockwidth + j] = 4139373.250000; float(dem_y(i, j));
h_dem_z[i * blockwidth + j] = 4393382.500000;float(dem_z(i, j));
h_demsloper_x[i * blockwidth + j] = 4393382.500000;float(demsloper_x(i, j));
h_demsloper_y[i * blockwidth + j] = 446.923950;float(demsloper_y(i, j));
h_demsloper_z[i * blockwidth + j] = -219.002213;float(demsloper_z(i, j));
h_dem_x[i * blockwidth + j] = -2028380.6250000; double(dem_x(i, j));
h_dem_y[i * blockwidth + j] = 4139373.250000; double(dem_y(i, j));
h_dem_z[i * blockwidth + j] = 4393382.500000;double(dem_z(i, j));
h_demsloper_x[i * blockwidth + j] = 4393382.500000;double(demsloper_x(i, j));
h_demsloper_y[i * blockwidth + j] = 446.923950;double(demsloper_y(i, j));
h_demsloper_z[i * blockwidth + j] = -219.002213;double(demsloper_z(i, j));
h_demcls[i * blockwidth + j] = clamap[80] ;// clamap[long(landcover(i, j))];
#else
h_dem_x[i * blockwidth + j] = float(dem_x(i, j));
h_dem_y[i * blockwidth + j] = float(dem_y(i, j));
h_dem_z[i * blockwidth + j] = float(dem_z(i, j));
h_demsloper_x[i * blockwidth + j] = float(demsloper_x(i, j));
h_demsloper_y[i * blockwidth + j] = float(demsloper_y(i, j));
h_demsloper_z[i * blockwidth + j] = float(demsloper_z(i, j));
h_dem_x[i * blockwidth + j] = double(dem_x(i, j));
h_dem_y[i * blockwidth + j] = double(dem_y(i, j));
h_dem_z[i * blockwidth + j] = double(dem_z(i, j));
h_demsloper_x[i * blockwidth + j] = double(demsloper_x(i, j));
h_demsloper_y[i * blockwidth + j] = double(demsloper_y(i, j));
h_demsloper_z[i * blockwidth + j] = double(demsloper_z(i, j));
h_demcls[i * blockwidth + j] = clamap[long(landcover(i, j))];
#endif
#endif
h_one[i * blockwidth + j] = 1;
}
}
HostToDevice((void*)h_dem_x, (void*)d_dem_x, sizeof(float) * newblokline * tempDemCols); // 复制 机器 -> GPU
HostToDevice((void*)h_dem_y, (void*)d_dem_y, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_dem_z, (void*)d_dem_z, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_x, (void*)d_demsloper_x, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_y, (void*)d_demsloper_y, sizeof(float) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_z, (void*)d_demsloper_z, sizeof(float) * newblokline * tempDemCols);
HostToDevice(h_one, d_one, sizeof(double) * newblokline * tempDemCols);
HostToDevice((void*)h_dem_x, (void*)d_dem_x, sizeof(double) * newblokline * tempDemCols); // 复制 机器 -> GPU
HostToDevice((void*)h_dem_y, (void*)d_dem_y, sizeof(double) * newblokline * tempDemCols);
HostToDevice((void*)h_dem_z, (void*)d_dem_z, sizeof(double) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_x, (void*)d_demsloper_x, sizeof(double) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_y, (void*)d_demsloper_y, sizeof(double) * newblokline * tempDemCols);
HostToDevice((void*)h_demsloper_z, (void*)d_demsloper_z, sizeof(double) * newblokline * tempDemCols);
HostToDevice((void*)h_demcls, (void*)d_demcls, sizeof(long) * newblokline * tempDemCols);
@ -688,6 +727,9 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
std::shared_ptr<double> h_temp_R(new double[PluseCount],delArrPtr);
#endif // __PRFDEBUG__
long pixelcount = newblokline * tempDemCols;
long startprfid = 0;
for (startprfid = 0; startprfid < pluseCount; startprfid = startprfid + echoblockline) {
@ -704,28 +746,28 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
for (long tempprfid = 0; tempprfid < templine; tempprfid++) {
{// 计算
long prfid = tempprfid + startprfid;
std::cout << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz").toStdString() << "]\t" << prfid<<"\t\t\t\t\t";
//std::cout << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz").toStdString() << "]\t" << prfid<<"\t start R\t\t\t\n";
// 天线位置
float antpx = sateOirbtNodes[prfid].Px;
float antpy = sateOirbtNodes[prfid].Py;
float antpz = sateOirbtNodes[prfid].Pz;
float antvx = sateOirbtNodes[prfid].Vx;
float antvy = sateOirbtNodes[prfid].Vy;
float antvz = sateOirbtNodes[prfid].Vz; //6
float antdirectx = sateOirbtNodes[prfid].AntDirecX;
float antdirecty = sateOirbtNodes[prfid].AntDirecY;
float antdirectz = sateOirbtNodes[prfid].AntDirecZ; // 9 天线指向
float antXaxisX = sateOirbtNodes[prfid].AntXaxisX;
float antXaxisY = sateOirbtNodes[prfid].AntXaxisY;
float antXaxisZ = sateOirbtNodes[prfid].AntXaxisZ;//12 天线坐标系
float antYaxisX = sateOirbtNodes[prfid].AntYaxisX;
float antYaxisY = sateOirbtNodes[prfid].AntYaxisY;
float antYaxisZ = sateOirbtNodes[prfid].AntYaxisZ;//15
float antZaxisX = sateOirbtNodes[prfid].AntZaxisX;
float antZaxisY = sateOirbtNodes[prfid].AntZaxisY;
float antZaxisZ = sateOirbtNodes[prfid].AntZaxisZ;//18
double antpx = sateOirbtNodes[prfid].Px;
double antpy = sateOirbtNodes[prfid].Py;
double antpz = sateOirbtNodes[prfid].Pz;
double antvx = sateOirbtNodes[prfid].Vx;
double antvy = sateOirbtNodes[prfid].Vy;
double antvz = sateOirbtNodes[prfid].Vz; //6
double antdirectx = sateOirbtNodes[prfid].AntDirecX;
double antdirecty = sateOirbtNodes[prfid].AntDirecY;
double antdirectz = sateOirbtNodes[prfid].AntDirecZ; // 9 天线指向
double antXaxisX = sateOirbtNodes[prfid].AntXaxisX;
double antXaxisY = sateOirbtNodes[prfid].AntXaxisY;
double antXaxisZ = sateOirbtNodes[prfid].AntXaxisZ;//12 天线坐标系
double antYaxisX = sateOirbtNodes[prfid].AntYaxisX;
double antYaxisY = sateOirbtNodes[prfid].AntYaxisY;
double antYaxisZ = sateOirbtNodes[prfid].AntYaxisZ;//15
double antZaxisX = sateOirbtNodes[prfid].AntZaxisX;
double antZaxisY = sateOirbtNodes[prfid].AntZaxisY;
double antZaxisZ = sateOirbtNodes[prfid].AntZaxisZ;//18
// 计算距离、局地入射角、增益
CUDARFPC_Caluation_R_Gain(
antpx, antpy, antpz, // 天线的坐标
@ -743,14 +785,18 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
NearRange, FarRange,
d_clsSigmaParam, clamapid,
d_factorj, PlusePoint,
d_R, // 输出距离
d_R, // 输出距离
//d_amp
d_PRFEcho_real, d_PRFEcho_imag,tempprfid // 输出振幅
);
//CUDA_PRF_CalFreqEcho(
// d_R, d_amp, pixelcount,
// d_factorj, PlusePoint,
// dx, NearRange,
// d_PRFEcho, tempprfid);
//std::cout << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz").toStdString() << "]\t" << prfid << "\t cal R\t\t\t\n";
//for (long fid = 0; fid < freqnum; fid++) {
// CUDA_PRF_GeneratorEcho(handle,d_R, d_amp, pixelcount,
// d_factorj, PlusePoint,
// NearRange, FarRange,
// d_PRFEcho_real, d_PRFEcho_imag, tempprfid);
//}
//std::cout << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz").toStdString() << "]\t" << prfid << "\t end Echo\t\t\t\n";
if (prfid % 100 == 0) {
std::cout << "[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz").toStdString() << "] dem:\t" << startline << "\t-\t" << startline + newblokline << "\t:\t pluse :\t" << prfid << " / " << pluseCount << std::endl;
@ -775,6 +821,10 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
#ifdef __PRFDEBUG__ && __PRFDEBUG_PRFINF__
testOutAmpArr("test_out_D_R.bin", h_temp_R.get(), pluseCount, 1);
#endif
#ifdef __PRFDEBUG__
break;
#endif // __PRFDEBUG__
@ -783,6 +833,8 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
std::cout << std::endl;
// 释放资源
cublasDestroy(handle);
// 地面数据释放
FreeCUDAHost(h_dem_x); FreeCUDADevice(d_dem_x);
@ -800,6 +852,8 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
FreeCUDAHost(h_factorj); FreeCUDADevice(d_factorj);
FreeCUDAHost(h_PRFEcho_real); FreeCUDADevice(d_PRFEcho_real);
FreeCUDAHost(h_PRFEcho_imag); FreeCUDADevice(d_PRFEcho_imag);
FreeCUDAHost(h_one); FreeCUDADevice(d_one);
FreeCUDAHost(h_temp); FreeCUDADevice(d_temp);
#endif

View File

@ -201,23 +201,28 @@ double AbstractSARSatelliteModel::getBandWidth()
QVector<double> AbstractSARSatelliteModel::getFreqList()
{
double bandwidth = this->getBandWidth(); // 频率范围
double centerFreq = this->getCenterFreq(); // 中心频率
double bandwidth = this->getBandWidth()/1e6; // ƵÂÊ·¶Î§
double centerFreq = this->getCenterFreq() / 1e6; // ÖÐÐÄÆµÂÊ
double nearR = this->getNearRange();
double farR = this->getFarRange();
// 计算分辨率
double Resolution = LIGHTSPEED / 2.0 / bandwidth; // 计算分辨率
long freqpoints = (farR - nearR) / Resolution + 1;
double minFreq = centerFreq - bandwidth / 2.0;// 最小频率
double maxFreq = minFreq + bandwidth;
double deltaFreq = bandwidth / (freqpoints - 1);
long long centerFreq_long = long long(centerFreq);
long long bandwidth_long = long long(bandwidth);
long long bandhalf = ceil(bandwidth_long / 2);
long long bw = bandhalf * 2;
double Resolution = LIGHTSPEED / 2.0/1e6 / bw;
long freqpoints = ceil((farR - nearR) / Resolution + 1);
Eigen::VectorXd freqs=Eigen::VectorXd::LinSpaced(freqpoints, -bandhalf, bandhalf);
freqs = freqs.array() + 1.0 * centerFreq_long;
QVector<double> freqlist(freqpoints);
for (long i = 0; i < freqpoints; i++) {
freqlist[i] = minFreq + i * deltaFreq;
freqlist[i] = freqs[i]*1e6;
}
return freqlist;
}