diff --git a/src/LAMPTool/BackScatterModel.cpp b/src/LAMPTool/BackScatterModel.cpp new file mode 100644 index 0000000..deab8bf --- /dev/null +++ b/src/LAMPTool/BackScatterModel.cpp @@ -0,0 +1,15 @@ +#include "LAMPToolAPI.h" +#include "BackScatterModel.h" +#include +#include + + +double LAMPTOOLAPI MuhlemanSimulationBackScatter(double incidentAngle) +{ + return 0.0133*cos(incidentAngle)/pow(sin(incidentAngle)+0.1*cos(incidentAngle), 3); +} + +Eigen::MatrixXd LAMPTOOLAPI MuhlemanSimulationBackScatter(Eigen::MatrixXd incidentAngle) +{ + return 0.0133 * (incidentAngle.array().cos()) / ((incidentAngle.array().sin()) + cos(incidentAngle.array().cos()*0.1).pow(3)); +} diff --git a/src/LAMPTool/BackScatterModel.h b/src/LAMPTool/BackScatterModel.h new file mode 100644 index 0000000..6fe85c3 --- /dev/null +++ b/src/LAMPTool/BackScatterModel.h @@ -0,0 +1,20 @@ +#pragma once +#ifndef _BACKSCATTERMODEL_H_ +#define _BACKSCATTERMODEL_H_ + + +#include "LAMPToolAPI.h" +#include +#include +#include "LAMPToolAPI.h" + +double LAMPTOOLAPI MuhlemanSimulationBackScatter(double incidentAngle); +Eigen::MatrixXd LAMPTOOLAPI MuhlemanSimulationBackScatter(Eigen::MatrixXd incidentAngle); + + +#endif + + + + + diff --git a/src/LAMPTool/BaseConstVariable.h b/src/LAMPTool/BaseConstVariable.h new file mode 100644 index 0000000..1aa9b4e --- /dev/null +++ b/src/LAMPTool/BaseConstVariable.h @@ -0,0 +1,38 @@ +#pragma once + +#ifndef BASECONSTVARIABLE_H +#define BASECONSTVARIABLE_H +#include "LAMPToolAPI.h" + +#include +#include + + +#define PI_180 180/3.141592653589793238462643383279 +#define T180_PI 3.141592653589793238462643383279/180 +#define LIGHTSPEED 299792458 + +#define Radians2Degrees(Radians) Radians*PI_180 +#define Degrees2Radians(Degrees) Degrees*T180_PI + + + +const double PI = 3.141592653589793238462643383279; +const double epsilon = 0.000000000000001; +const double pi = 3.14159265358979323846; +const double d2r = pi / 180; +const double r2d = 180 / pi; + +const double a = 6378137.0; //椭球长半轴 +const double ae = 6378137.0; //椭球长半轴 +const double ee = 0.0818191910428;// 第一偏心率 +const double f_inverse = 298.257223563; //扁率倒数 +const double b = a - a / f_inverse; +const double eSquare = (a * a - b * b) / (a * a); +const double e = sqrt(eSquare); +const double earth_Re = 6378136.49; +const double earth_Rp = (1 - 1 / f_inverse) * earth_Re; +const double earth_We = 0.000072292115; + + +#endif \ No newline at end of file diff --git a/src/LAMPTool/BaseTool.cpp b/src/LAMPTool/BaseTool.cpp new file mode 100644 index 0000000..05045a1 --- /dev/null +++ b/src/LAMPTool/BaseTool.cpp @@ -0,0 +1,256 @@ +#pragma once +#include "BaseTool.h" +/// +/// +//#define EIGEN_USE_MKL_ALL +//#define EIGEN_VECTORIZE_SSE4_2 +//#include + +#include "referenceHeader.h" +#include +#include +#include +#include +//#include +#include +#include + +#include < io.h > +#include < stdio.h > +#include < stdlib.h > +#include +#include +#include +//#include +#include //#include "ogrsf_frmts.h" + +#include +#include +#include "GeoOperator.h" + + +#include "baseTool.h" +using namespace std; +using namespace Eigen; + + + + +QString LAMPTOOLAPI getCurrentTimeString() { + struct tm ConversionTime; + std::time_t t = std::time(NULL); + char mbstr[100]; + _localtime64_s(&ConversionTime, &t); + std::strftime(mbstr, sizeof(mbstr), "%Y-%m-%d %H:%M:%S", &ConversionTime); + QString strTime = mbstr; + return strTime; +} + +QString LAMPTOOLAPI getCurrentShortTimeString() { + struct tm ConversionTime; + std::time_t t = std::time(NULL); + char mbstr[100]; + _localtime64_s(&ConversionTime, &t); + std::strftime(mbstr, sizeof(mbstr), "%Y-%m-%d %H:%M:%S", &ConversionTime); + QString strTime = mbstr; + return strTime; +} + +std::vector LAMPTOOLAPI splitString(const QString& str, char delimiter) +{ + QStringList tokens = str.split(delimiter); + return convertQStringListToStdVector(tokens); +} + + +complex LAMPTOOLAPI Cubic_Convolution_interpolation(double u, double v, Eigen::MatrixX> img) +{ + if (img.rows() != 4 || img.cols() != 4) { + throw exception("the size of img's block is not right"); + } + // 斤拷锟斤拷模锟斤拷 + Eigen::MatrixX> wrc(1, 4);// 使锟斤拷 complex 斤拷锟斤拷要原斤拷为锟剿伙拷取值 + Eigen::MatrixX> wcr(4, 1);// + for (int i = 0; i < 4; i++) { + wrc(0, i) = Cubic_kernel_weight(u + 1 - i); // u+1,u,u-1,u-2 + wcr(i, 0) = Cubic_kernel_weight(v + 1 - i); + } + + Eigen::MatrixX> interValue = wrc * img * wcr; + return interValue(0, 0); +} + +complex LAMPTOOLAPI Cubic_kernel_weight(double s) +{ + s = abs(s); + if (s <= 1) { + return complex(1.5 * s * s * s - 2.5 * s * s + 1, 0); + } + else if (s <= 2) { + return complex(-0.5 * s * s * s + 2.5 * s * s - 4 * s + 2, 0); + } + else { + return complex(0, 0); + } +} + +/// +/// p11 p12 -- x +/// p0(u,v) +/// p21 p22 +/// | +/// y +/// p11(0,0) +/// p21(0,1) +/// P12(1,0) +/// p22(1,1) +/// +/// x,y,z +/// x,y,z +/// x,y,z +/// x,y,z +/// x,y,z +/// +double LAMPTOOLAPI Bilinear_interpolation(Landpoint p0, Landpoint p11, Landpoint p21, Landpoint p12, Landpoint p22) +{ + + return p11.ati * (1 - p0.lon) * (1 - p0.lat) + + p12.ati * p0.lon * (1 - p0.lat) + + p21.ati * (1 - p0.lon) * p0.lat + + p22.ati * p0.lon * p0.lat; +} + + + +bool LAMPTOOLAPI onSegment(Point_3d Pi, Point_3d Pj, Point_3d Q) +{ + if ((Q.x - Pi.x) * (Pj.y - Pi.y) == (Pj.x - Pi.x) * (Q.y - Pi.y) //叉乘 + //保证Q点坐标在pi,pj之间 + && min(Pi.x, Pj.x) <= Q.x && Q.x <= max(Pi.x, Pj.x) + && min(Pi.y, Pj.y) <= Q.y && Q.y <= max(Pi.y, Pj.y)) + return true; + else + return false; +} + +Point_3d LAMPTOOLAPI invBilinear(Point_3d p, Point_3d a, Point_3d b, Point_3d c, Point_3d d) +{ + Point_3d res; + + Point_3d e = b - a; + Point_3d f = d - a; + Point_3d g = a - b + c - d; + Point_3d h = p - a; + + double k2 = cross2d(g, f); + double k1 = cross2d(e, f) + cross2d(h, g); + double k0 = cross2d(h, e); + double u, v; + // if edges are parallel, this is a linear equation + if (abs(k2) < 0.001) + { + v = -k0 / k1; + u = (h.x - f.x * v) / (e.x + g.x * v); + p.z = a.z + (b.z - a.z) * u + (d.z - a.z) * v + (a.z - b.z + c.z - d.z) * u * v; + return p; + } + // otherwise, it's a quadratic + else + { + float w = k1 * k1 - 4.0 * k0 * k2; + if (w < 0.0) { + // 可能在边界上 + if (onSegment(a, b, p)) { + Point_3d tt = b - a; + Point_3d ttpa = p - a; + double scater = ttpa / tt; + if (scater < 0 || scater>1) { return { -9999,-9999,-9999 }; } + p.z = a.z + scater * tt.z; + return p; + } + else if (onSegment(b, c, p)) { + Point_3d tt = c - b; + Point_3d ttpa = p - b; + double scater = ttpa / tt; + if (scater < 0 || scater>1) { return { -9999,-9999,-9999 }; } + p.z = b.z + scater * tt.z; + return p; + } + else if (onSegment(c, d, p)) { + Point_3d tt = d - c; + Point_3d ttpa = p - c; + double scater = ttpa / tt; + if (scater < 0 || scater>1) { return { -9999,-9999,-9999 }; } + p.z = c.z + scater * tt.z; + return p; + } + else if (onSegment(d, a, p)) { + Point_3d tt = a - d; + Point_3d ttpa = p - d; + double scater = ttpa / tt; + if (scater < 0 || scater>1) { return { -9999,-9999,-9999 }; } + p.z = d.z + scater * tt.z; + return p; + } + + return { -9999,-9999,-9999 }; + } + else { + w = sqrt(w); + + float ik2 = 0.5 / k2; + float v = (-k1 - w) * ik2; + float u = (h.x - f.x * v) / (e.x + g.x * v); + + if (u < 0.0 || u>1.0 || v < 0.0 || v>1.0) + { + v = (-k1 + w) * ik2; + u = (h.x - f.x * v) / (e.x + g.x * v); + } + p.z = a.z + (b.z - a.z) * u + (d.z - a.z) * v + (a.z - b.z + c.z - d.z) * u * v; + return p; + } + } + p.z = a.z + (b.z - a.z) * u + (d.z - a.z) * v + (a.z - b.z + c.z - d.z) * u * v; + return p; +} + +double LAMPTOOLAPI sind(double degree) +{ + return sin(degree * d2r); +} + +double LAMPTOOLAPI cosd(double d) +{ + return cos(d * d2r); +} + + +string LAMPTOOLAPI Convert(float Num) +{ + ostringstream oss; + oss << Num; + string str(oss.str()); + return str; +} + +QString LAMPTOOLAPI JoinPath(const QString& path, const QString& filename) +{ + QDir dir(path); + + // Ensure the path ends with the appropriate separator + if (!QDir::isAbsolutePath(path)) + dir.makeAbsolute(); + + return dir.filePath(filename); +} +std::vector LAMPTOOLAPI convertQStringListToStdVector(const QStringList& qStringList) +{ + std::vector stdVector; + + for (const QString& str : qStringList) { + stdVector.push_back(str); + } + + return stdVector; +}; \ No newline at end of file diff --git a/src/LAMPTool/BaseTool.h b/src/LAMPTool/BaseTool.h new file mode 100644 index 0000000..54d76a0 --- /dev/null +++ b/src/LAMPTool/BaseTool.h @@ -0,0 +1,94 @@ +#pragma once +#pragma once +#ifndef BASETOOL_H +#define BASETOOL_H + +/// +/// 基本类、基本函数 +/// +//#define EIGEN_USE_MKL_ALL +//#define EIGEN_VECTORIZE_SSE4_2 +//#include + +//#include +#include +#include +#include +#include +#include +#include +#include +#include "referenceHeader.h" +#include +#include +#include +#include +#include +#include "GeoOperator.h" +#include +#include + + +using namespace std; +using namespace Eigen; + + +///////////////////////////////////// 运行时间打印 ////////////////////////////////////////////////////////// + + +QString LAMPTOOLAPI getCurrentTimeString(); +QString LAMPTOOLAPI getCurrentShortTimeString(); + +std::vector LAMPTOOLAPI splitString(const QString& str, char delimiter); + + + +/////////////////////////////// 基本图像类 结束 ////////////////////////////////////////////////////////// + +string LAMPTOOLAPI Convert(float Num); +QString LAMPTOOLAPI JoinPath(const QString& path, const QString& filename); + +////////////////////////////// 坐标部分基本方法 ////////////////////////////////////////// + + +////////////////////////////// 坐标部分基本方法 ////////////////////////////////////////// + + +////////////////////////////// 插值 //////////////////////////////////////////// + +complex LAMPTOOLAPI Cubic_Convolution_interpolation(double u, double v, Eigen::MatrixX> img); + +complex LAMPTOOLAPI Cubic_kernel_weight(double s); + +double LAMPTOOLAPI Bilinear_interpolation(Landpoint p0, Landpoint p11, Landpoint p21, Landpoint p12, Landpoint p22); + +bool LAMPTOOLAPI onSegment(Point_3d Pi, Point_3d Pj, Point_3d Q); + +Point_3d LAMPTOOLAPI invBilinear(Point_3d p, Point_3d a, Point_3d b, Point_3d c, Point_3d d); + + + +// +// WGS84 到J2000 坐标系的变换 +// 参考网址：https://blog.csdn.net/hit5067/article/details/116894616 +// 资料网址：http://celestrak.org/spacedata/ +// 参数文件： +// a. Earth Orientation Parameter 文件： http://celestrak.org/spacedata/EOP-Last5Years.csv +// b. Space Weather Data 文件： http://celestrak.org/spacedata/SW-Last5Years.csv +// 备注：上述文件是自2017年-五年内 +/** +在wgs84 坐标系转到J2000 坐标系 主要 涉及到坐标的相互转换。一般给定的WGS坐标为 给定时刻的 t ,BLH +转换步骤： +step 1: WGS 84 转换到协议地球坐标系 +step 2: 协议地球坐标系 转换为瞬时地球坐标系 +step 3: 瞬时地球坐标系 转换为 瞬时真天球坐标系 +step 4: 瞬时真天球坐标系 转到瞬时平天球 坐标系 +step 5: 瞬时平天球坐标系转换为协议天球坐标系（J2000） +**/ + + +double LAMPTOOLAPI sind(double degree); + +double LAMPTOOLAPI cosd(double d); + +#endif \ No newline at end of file diff --git a/src/LAMPTool/FEKOBaseToolClass.cpp b/src/LAMPTool/FEKOBaseToolClass.cpp new file mode 100644 index 0000000..934fe3c --- /dev/null +++ b/src/LAMPTool/FEKOBaseToolClass.cpp @@ -0,0 +1,1040 @@ +#include "FEKOBaseToolClass.h" +#include +#include "BaseConstVariable.h" +#include +#include +#include +#include "ImageOperatorBase.h" +#include "SARImageBase.h" +#include "SARBaseTool.h" +#include "FileOperator.h" +#include + + + + + + + + + +QString LAMPTOOLAPI FEKOBase::FEKOImageModeenumToString(FEKOBase::FEKOImageMode mode) +{ + switch (mode) { + case Strip: + return "STRIP"; + case Scane: + return "SCANE"; + case ISAR: + return "ISAR"; + case CircleSAR: + return "CircleSAR"; + } + // 如果枚举值不匹配任何字符串，则返回一个默认值，或者抛出异常 + return "UNKNOW"; // 默认返回Unknown +} + +FEKOBase::FEKOCoordinateSystem LAMPTOOLAPI FEKOBase::FEKOCoordinateSystemString2Enum(QString str) +{ + if (str.toUpper() == "SPHERICAL") { return FEKOBase::Spherical; } + else if (str.toUpper() == "CARTESIAN") { return FEKOBase::Cartesian; } + else { + return FEKOBase::UNKONWFEKOCOORDINATESYSTEM; + } + +} + +QString LAMPTOOLAPI FEKOBase::QString2FEKOCoordinateSystem(FEKOBase::FEKOCoordinateSystem mode) +{ + switch (mode) { + case FEKOBase::Spherical: + return "SPHERICAL"; + case FEKOBase::Cartesian: + return "CARTESIAN"; + default: + return "UNKONWFEKOCOORDINATESYSTEM"; + + } +} + +FEKOBase::FEKOImageMode LAMPTOOLAPI FEKOBase::FEKOImageModeString2Enum(QString str) +{ + if (str.toUpper() == "STRIP") { + return FEKOBase::FEKOImageMode::Strip; + } + else if (str.toUpper() == "SCANE") { + return FEKOBase::FEKOImageMode::Scane; + } + else if (str.toUpper() == "ISAR") { + return FEKOBase::FEKOImageMode::ISAR; + } + else if (str.toUpper() == "CIRCLESAR") { + return FEKOBase::FEKOImageMode::CircleSAR; + } + // 如果输入的字符串不匹配任何枚举值，则返回一个默认值，或者抛出异常 + return FEKOBase::FEKOImageMode::UNKNOW; // 默认返回UNKNOW +} + +FEKOBase::freqParams LAMPTOOLAPI FEKOBase::getFreqSetting(double centerFreq, double resolution, double bandWidth, double scenceRange, bool isResolution) { + + FEKOBase::freqParams result{ 0,0,0 }; + { + if (isResolution) { + bandWidth = 0.299792458 / 2 / resolution; // 计算带宽 + } + else { + resolution = 0.299792458 / 2 / bandWidth; // 计算分辨率 + } + } + size_t samplePoint = std::ceil((scenceRange) / resolution) + 1; + result.startfreqs = centerFreq - 0.5 * bandWidth; + result.endfreqs = centerFreq + 0.5 * bandWidth; + result.freqpoint = samplePoint; + return result; +} + +FEKOBase::FEKOSatelliteParams LAMPTOOLAPI FEKOBase::createFEKOSatelliteParams(double Px, double Py, double Pz, double Vx, double Vy, double Vz, double incidenceAngle, double AzAngle, double theta, double phi, bool isRight, size_t PRFIdx) +{ + FEKOBase::FEKOSatelliteParams result{ + PRFIdx, + FEKOBase::SatelliteState{ + FEKOBase::SatellitePosition{Px,Py,Pz}, + FEKOBase::SatelliteVelocity{Vx,Vy,Vz} + }, + incidenceAngle,AzAngle, + FEKOBase::FEKOantPitionDirect{Px,Py,Pz,theta,phi}, + isRight + }; + return result; +} + +FEKOBase::FEKOSatelliteParams LAMPTOOLAPI FEKOBase::createFEKOSatelliteParams(SatelliteState pose, double incidenceAngle, double AzAngle, FEKOantPitionDirect antpos, size_t PRFIdx) +{ + FEKOBase::FEKOSatelliteParams result{ + PRFIdx,pose,incidenceAngle,AzAngle,antpos + }; + return result; +} + +FEKOBase::SatelliteState LAMPTOOLAPI FEKOBase::FEKOSatelliteParams2SatelliteState(FEKOSatelliteParams parmas) +{ + return parmas.pose; +} + +FEKOBase::FEKOantPitionDirect LAMPTOOLAPI FEKOBase::FEKOSatelliteParams2FEKOantPitionDirect(FEKOSatelliteParams parmas) +{ + return parmas.antpos; +} + +TopoDS_Shape LAMPTOOLAPI FEKOBase::SatellitePos2FEKOAntPos(SatelliteState satepos, double incidenceAngle, double AzAngle, bool isRIGHT, FEKOantPitionDirect* antposition_Direct, TopoDS_Shape antModel) +{ + incidenceAngle = incidenceAngle * M_PI / 180; // 转换为弧度 + AzAngle = AzAngle * M_PI / 180; // 扫描角度转换为弧度 + TopoDS_Shape tempShape; + if (!antModel.IsNull()) { + // 创建一个复制 + BRepBuilderAPI_Copy copyBuilder(antModel); + tempShape = copyBuilder.Shape(); + } + + double Sx = satepos.pos.Px; + double Sy = satepos.pos.Py; + double Sz = satepos.pos.Pz; + double Vx = satepos.vel.Vx; + double Vy = satepos.vel.Vy; + double Vz = satepos.vel.Vz; + //qDebug() << u8"// 0. 卫星状态矢量"; + gp_Vec SateVelocity(gp_Pnt(0, 0, 0), gp_Pnt(Vx, Vy, Vz)); // 卫星速度矢量 + gp_Vec SatePosition(gp_Pnt(0, 0, 0), gp_Pnt(Sx, Sy, Sz)); // 卫星位置矢量 + gp_Vec sch_X(gp_Pnt(0, 0, 0), gp_Pnt(1.0, 0.0, 0.0)); + gp_Vec sch_Y(gp_Pnt(0, 0, 0), gp_Pnt(0.0, 1.0, 0.0)); // 飞行方向 y + gp_Vec sch_Z(gp_Pnt(0, 0, 0), gp_Pnt(0.0, 0.0, 1.0)); // 雷达照射方向 + +// qDebug() << u8"// 1. 卫星坐标系 sch -->飞行坐标系 fly 保证Z轴指向不变化"; + gp_Vec sch_X_fly; + gp_Vec sch_Y_fly; // 飞行向 + gp_Vec sch_Z_fly; + + { + CartesianCoordinates Sp_car{ sch_Y.X(),sch_Y.Y(), sch_Y.Z() }; + CartesianCoordinates tp_car{ SateVelocity.X(),SateVelocity.Y(), SateVelocity.Z() }; + + SphericalCoordinates sp_sph = cartesianToSpherical(Sp_car); + SphericalCoordinates tp_sph = cartesianToSpherical(tp_car); + + // 计算 + double delta_theta = tp_sph.theta - sp_sph.theta; + double delta_phi = tp_sph.phi - sp_sph.phi; + + gp_Trsf rotationTransform_theta; + rotationTransform_theta.SetRotation(gp_Ax1(gp_Pnt(0, 0, 0), gp_Dir(0, 1, 0)), delta_theta); // 绕 y 轴旋转 + sch_X_fly = sch_X.Transformed(rotationTransform_theta); + sch_Y_fly = sch_Y.Transformed(rotationTransform_theta); + sch_Z_fly = sch_Z.Transformed(rotationTransform_theta); + + gp_Trsf rotationTransform_phi; + rotationTransform_phi.SetRotation(gp_Ax1(gp_Pnt(0, 0, 0), gp_Dir(0, 0, 1)), delta_phi); // 绕 z 轴旋转 + sch_X_fly = sch_X_fly.Transformed(rotationTransform_phi); + sch_Y_fly = sch_Y_fly.Transformed(rotationTransform_phi); + sch_Z_fly = sch_Z_fly.Transformed(rotationTransform_phi); + + if (!tempShape.IsNull()) { + BRepBuilderAPI_Transform shapeTransform_theta(tempShape, rotationTransform_theta); + tempShape = shapeTransform_theta.Shape(); + BRepBuilderAPI_Transform shapeTransform_phi(tempShape, rotationTransform_phi); + tempShape = shapeTransform_phi.Shape(); + } + } + +// qDebug() << u8"// 2. 根据方位角与俯仰角，调整雷达姿态 --- 此步确定 Y 轴已经调整完毕 "; +// qDebug() << u8"// 2. a 根据 incidenceAngle 调整 Z 的指向，需要先根据左右视判断出雷达照射方向"; + + { // 按照 Y 轴进行计算 + incidenceAngle = M_PI - incidenceAngle; // 下Z轴 + Standard_Real ZRotation_angle = 0; + if (isRIGHT) { // 右视 + + ZRotation_angle = incidenceAngle; + } + else { // 左视 - + ZRotation_angle = -1 * incidenceAngle; + } + gp_Trsf fly_incidence_trsf; + fly_incidence_trsf.SetRotation(gp_Ax1(gp_Pnt(0, 0, 0), sch_Y_fly), ZRotation_angle); // 计算入射角 + + sch_X_fly = sch_X_fly.Transformed(fly_incidence_trsf); + sch_Y_fly = sch_Y_fly.Transformed(fly_incidence_trsf); + sch_Z_fly = sch_Z_fly.Transformed(fly_incidence_trsf); // 此时获取雷达的照射方向 + + if (!tempShape.IsNull()) { + BRepBuilderAPI_Transform shapeTransform_incidence(tempShape, fly_incidence_trsf); // fly_Z --> incidenceAngle + tempShape = shapeTransform_incidence.Shape(); + } + + } + +// qDebug() << u8"// 2.b 根据 AzAngle 调整 Z 的指向 "; +// qDebug() << u8"// Az 应该变化在 YoZ， 也就是绕 X 轴 旋转；"; + + { // 沿Z 轴进行计算 + gp_Trsf fly_AZ_trsf; + fly_AZ_trsf.SetRotation(gp_Ax1(gp_Pnt(0, 0, 0), sch_X_fly), AzAngle); // 计算方位角变化 + + sch_X_fly = sch_X_fly.Transformed(fly_AZ_trsf); + sch_Y_fly = sch_Y_fly.Transformed(fly_AZ_trsf); + sch_Z_fly = sch_Z_fly.Transformed(fly_AZ_trsf); // 此时获取雷达的照射方向 + + if (!tempShape.IsNull()) { + BRepBuilderAPI_Transform shapeTransform_AZ(tempShape, fly_AZ_trsf); // fly_Z --> incidenceAngle + tempShape = shapeTransform_AZ.Shape(); // 旋转卫星模型 + } + + } +// qDebug() << u8"// 2.c 参数记录 "; + { + CartesianCoordinates Sp_car{ 0,0,1 }; + CartesianCoordinates tp_car{ sch_Z_fly.X(),sch_Z_fly.Y(), sch_Z_fly.Z() }; + + SphericalCoordinates sp_sph = cartesianToSpherical(Sp_car); + SphericalCoordinates tp_sph = cartesianToSpherical(tp_car); // 计算phi角度 + + double delta_theta = tp_sph.theta - sp_sph.theta; + double delta_phi = tp_sph.phi - sp_sph.phi; + + antposition_Direct->x = Sx; //-- FEKO 设置 + antposition_Direct->y = Sy; //-- FEKO 设置 + antposition_Direct->z = Sz; //-- FEKO 设置 + antposition_Direct->theta = delta_theta / M_PI * 180; //-- FEKO 设置 theta + antposition_Direct->phi = delta_phi / M_PI * 180;//--- FEKO 设置 phi + + } +// qDebug() << u8"//3. 执行模型平移"; + { + if (!tempShape.IsNull()) { + gp_Trsf fly_move_trsf; + fly_move_trsf.SetTranslation(SatePosition); + BRepBuilderAPI_Transform shapeTransform_move(tempShape, fly_move_trsf); // fly_Z --> incidenceAngle + tempShape = shapeTransform_move.Shape(); // 平移卫星模型 + + } + } + + qDebug() << QString(u8"雷达坐标与姿态（X,Y,Z,theta,phi): （%1,%2,%3,%4,%5)").arg(antposition_Direct->x) + .arg(antposition_Direct->y) + .arg(antposition_Direct->z) + .arg(antposition_Direct->theta) + .arg(antposition_Direct->phi); + + return tempShape; +} + +bool LAMPTOOLAPI FEKOBase::compareElectricFieldDataInFreq(const ElectricFieldData& a, const ElectricFieldData& b) +{ + return a.frequency < b.frequency; +} + +bool LAMPTOOLAPI FEKOBase::comparePRFPluseDataInPRFIdx(const PRFPluseData& a, const PRFPluseData& b) +{ + return a.prfidx < b.prfidx; +} + +LAMPTOOLAPI FEKOBase::NearFieldEchoCSVParser::NearFieldEchoCSVParser() +{ +} + +LAMPTOOLAPI FEKOBase::NearFieldEchoCSVParser::~NearFieldEchoCSVParser() +{ +} + +bool LAMPTOOLAPI FEKOBase::NearFieldEchoCSVParser::checkPRFModel() +{ + + qDebug() << u8"正在检查是否可以采用脉冲计数模式 configuration Name ：PRF_{脉冲计数}"; + + QRegExp regex("(\\d+)"); + // 判断是否需要按照脉冲次序排序 + bool usePRFbool = true; + QMap configName_prfidx; + // 使用迭代器遍历 QMap + QMap>::const_iterator it; + for (it = this->electricFieldDataList.constBegin(); it != this->electricFieldDataList.constEnd(); ++it) { + for (size_t i = 0; i < this->electricFieldDataList[it.key()].size(); i++) { + QString value = this->electricFieldDataList[it.key()][i].configurationName; + if (value.indexOf("PRF") != -1) { + if (regex.indexIn(value) != -1) { + // 提取匹配到的部分 + QString matchedText = regex.cap(1); + + size_t PRFidx = matchedText.toInt(&usePRFbool); + if (usePRFbool) { + configName_prfidx.insert(this->electricFieldDataList[it.key()][i].configurationName, PRFidx); + } + else { + qDebug() << u8"不能采用脉冲计数模型"; + QMessageBox::information(nullptr, u8"警告", u8"无法根据Configuration Name识别脉冲顺序，按照读取顺序，计数PRF", QMessageBox::Ok); + this->usePRFCountMode = false; + return false; + } + } + } + else { + qDebug() << u8"不能采用脉冲计数模型"; + QMessageBox::information(nullptr, u8"警告", u8"无法根据Configuration Name识别脉冲顺序，按照读取顺序，计数PRF", QMessageBox::Ok); + this->usePRFCountMode = false; + return false; + } + } + } + + // 更新prfidx + for (it = this->electricFieldDataList.constBegin(); it != this->electricFieldDataList.constEnd(); ++it) { + if (!(this->electricFieldDataList[it.key()].size() > 0)) { + QMessageBox::warning(nullptr, u8"错误", u8"发现空脉冲：" + it.key()); + qWarning() << u8"发现空脉冲：" + it.key(); + return false; + } + else {} + + QString value = this->electricFieldDataList[it.key()][0].configurationName; + + for (size_t i = 0; i < this->electricFieldDataList[it.key()].size(); i++) { + if (value != this->electricFieldDataList[it.key()][0].configurationName) { + QMessageBox::warning(nullptr, u8"错误", u8"脉冲解析错误，名称不一致：" + it.key() + " | " + value); + qWarning() << u8"脉冲解析错误，名称不一致：" + it.key() + " | " + value; + } + else {} + size_t tempprfidx = configName_prfidx[value]; + this->electricFieldDataList[it.key()][i].prfidx = tempprfidx; + } + } + + qDebug() << u8"结束检查是否可以采用脉冲计数模式 configuration Name ：PRF_{脉冲计数}"; + return true; +} + +bool LAMPTOOLAPI FEKOBase::NearFieldEchoCSVParser::resizePRFPluse() +{ + qDebug() << u8"根据数据文件，重新整理成脉冲形式"; + + QMap prfPluseMap; // 采用字典映射方式 + + qDebug() << u8"QMap> ==> QMap"; + double px; // 2. 计算姿态，近场点坐标 + double py; + double pz; + double R; + double theta; // x0y + double phi; // y0z + size_t PRFidx; + + for (QMap>::const_iterator it = this->electricFieldDataList.constBegin(); it != this->electricFieldDataList.constEnd(); ++it) { + PRFPluseData temp; + std::vector tempData = this->electricFieldDataList[it.key()]; // 脉冲已经排序号 + std::sort(tempData.begin(), tempData.end(), compareElectricFieldDataInFreq); + temp.electricFieldDataList = tempData; // 按照脉冲排序 + + temp.freqstart = tempData[0].frequency; // 1. 频率 + temp.freqend = tempData[tempData.size() - 1].frequency; + temp.freqpoints = tempData.size(); + + px = tempData[0].origin.x; // 2. 近场接收点计算姿态，近场点坐标 + py = tempData[0].origin.y; + pz = tempData[0].origin.z; + + + R = tempData[0].radius; + theta = tempData[0].theta; // x0y + phi = tempData[0].phi; // y0z + + + + + PRFidx = tempData[0].prfidx; + qDebug() << u8"分析姿态:" + it.key(); + for (size_t i = 0; i < tempData.size(); i++) { + if (abs(px - tempData[i].origin.x) > 0.001 || abs(py - tempData[i].origin.y) > 0.001 || + abs(pz - tempData[i].origin.z) > 0.001 || abs(R - tempData[i].radius) > 0.001 || + abs(theta - tempData[i].theta) > 0.001 || abs(phi - tempData[i].phi) > 0.001) + { + qDebug() << u8"发现 configuration Name:" + it.key() + u8" , 发现问题，请检查"; + QMessageBox::warning(nullptr, u8"警告", u8"发现 configuration Name:" + it.key() + u8" , 发现问题，请检查"); + return false; + } + else {} + } + // 计算增量 R, theta, phi ==> delta_R,delta_theta,delta_phi + theta = theta * M_PI / 180; + phi = phi * M_PI / 180; + + px = px + R * std::sin(theta) * std::cos(phi); + py = py + R * std::sin(theta) * std::sin(phi); + pz = pz + R * std::cos(theta); + + temp.px = px; + temp.py = py; + temp.pz = pz; + temp.prfidx = PRFidx; + + + prfPluseMap.insert(it.key(), temp); + } + + + qDebug() << u8"判断是所有回波均一性"; + this->prfData = std::vector(0); + size_t freqpointnum = prfPluseMap.first().freqpoints; + for (QMap::const_iterator it = prfPluseMap.constBegin(); it != prfPluseMap.constEnd(); ++it) { + if (it.value().freqpoints != freqpointnum) { + qDebug() << u8"回波数据频点数不一致，无法进行脉冲计数"; + QMessageBox::information(nullptr, u8"警告", u8"回波数据频点数不一致，无法进行脉冲计数", QMessageBox::Ok); + this->usePRFCountMode = false; + return false; + } + this->prfData.push_back(it.value()); + } + this->freqPoints = freqpointnum; + std::sort(this->prfData.begin(), this->prfData.end(), comparePRFPluseDataInPRFIdx); + + this->freqStart = this->prfData[0].freqstart; + this->freqEnd = this->prfData[0].freqend; + qDebug() << u8"根据数据文件，结束整理成脉冲形式"; + return true; + +} + +bool LAMPTOOLAPI FEKOBase::NearFieldEchoCSVParser::parseCSV(const QString& filePath) +{ + + std::vector dataList; + + QFile file(filePath); + if (file.open(QIODevice::ReadOnly | QIODevice::Text)) { + QTextStream in(&file); + + // 读取标题行 + QString headerLine = in.readLine(); + QStringList headers = headerLine.split(QRegularExpression(",(?!\\s)")); + + while (!in.atEnd()) { + QString line = in.readLine(); + + QStringList fields = line.split(QRegularExpression(",(?!\\s)")); + + ElectricFieldData data; + // 解析数据根据列名 + for (int i = 0; i < headers.size(); ++i) { + QString header = headers[i].trimmed(); + QString value = fields[i];// .trimmed(); + + if (header == "File Type") { + data.fileType = value; + } + else if (header == "File Format") { + data.fileFormat = value; + } + else if (header == "Source") { + data.source = value; + } + else if (header == "Date") { + data.date = value; + } + else if (header == "Radius") { // 半径 + data.radius = value.toDouble(); + } + else if (header == "Theta") { // theta 角 + data.theta = value.toDouble(); + } + else if (header == "Phi") {// phi 角 + data.phi = value.toDouble(); + } + else if (header == "Re(Er)") { + data.reEr = value.toDouble(); + } + else if (header == "Im(Er)") { + data.imEr = value.toDouble(); + } + else if (header == "Re(Etheta)") { + data.reEtheta = value.toDouble(); + } + else if (header == "Im(Etheta)") { + data.imEtheta = value.toDouble(); + } + else if (header == "Re(Ephi)") { + data.reEphi = value.toDouble(); + } + else if (header == "Im(Ephi)") { + data.imEphi = value.toDouble(); + } + else if (header == "Configuration Name") { // 脉冲次序名称 + data.configurationName = value.trimmed(); + } + else if (header == "Request Name") { + data.requestName = value.trimmed(); + } + else if (header == "Frequency") { + data.frequency = value.toDouble(); + } + else if (header == "Coordinate System") { + data.coordinateSystem = value; + } + else if (header == "Origin") { // 坐标参考原点 + QStringList vls = value.replace("\"", "").replace("(", "").replace(")", "").split(","); + if (vls.length() < 3) { + qDebug() << "******************** ERROR INFO ******************************************"; + qDebug() << "origin point has error ," << value; + qDebug() << "line: " << line; + qDebug() << "head: " << header; + qDebug() << "***************************************************************************"; + QMessageBox::information(nullptr, u8"警告", u8"无法解析Origin 参数，请检查", QMessageBox::Ok); + return false; + } + data.origin = Point_3d{ vls.at(0).toDouble(), + vls.at(1).toDouble(), + vls.at(2).toDouble() }; + } + else if (header == "Num Radius Samples") { + data.numRadiusSamples = value.toInt(); + } + else if (header == "Num Theta Samples") { + data.numThetaSamples = value.toInt(); + } + else if (header == "Num Phi Samples") { + data.numPhiSamples = value.toInt(); + } + else if (header == "Result Type") { + data.resultType = value; + } + else if (header == "Num Header Lines") { + data.numHeaderLines = value.toInt(); + } + // 添加其他字段的解析... + } + dataList.push_back(data); + } + + file.close(); + } + else { + qDebug() << "Failed to open the file:" << file.errorString(); + return false; + } + + //this->electricFieldDataList = dataList; + + qDebug() << u8"文件读取结果，对文件结果进行重新整理与检查"; + //获取脉冲最小值，最大值，如果是 使用软件生产的 configname 应该为 PRF_1 PRF_2 这里通过判断 1,2 来确定脉冲序号 + QSet configNameSet; + size_t curPRF = 0; + for (int i = 0; i < dataList.size(); i++) { + if (!configNameSet.contains(dataList[i].configurationName)) { + curPRF = curPRF + 1; + } + dataList[i].prfidx = curPRF; + configNameSet.insert(dataList[i].configurationName); + } + + // 初始化 + for (const QString& elementconfigName : configNameSet) { + this->electricFieldDataList.insert(elementconfigName, std::vector(0)); + } + + // 数据初始化 + for (size_t i = 0; i < dataList.size(); i++) { + this->electricFieldDataList[dataList[i].configurationName].push_back(dataList[i]); + } + // 每个Pluse 进行排序,按照频率： 小 -> 大， + QMap>::const_iterator it; + for (it = this->electricFieldDataList.constBegin(); it != this->electricFieldDataList.constEnd(); ++it) { + std::sort(this->electricFieldDataList[it.key()].begin(), this->electricFieldDataList[it.key()].end(), compareElectricFieldDataInFreq); + } + + // 检查每个electricFieldDataList 是否按照频率进行排序 + for (it = this->electricFieldDataList.constBegin(); it != this->electricFieldDataList.constEnd(); ++it) { + for (size_t i = 0; i < this->electricFieldDataList[it.key()].size() - 1; i++) { + if (this->electricFieldDataList[it.key()][i].frequency > this->electricFieldDataList[it.key()][i + 1].frequency) { + qDebug() << u8"发现 configuration Name:" + it.key() + u8" , 发现问题，请检查"; + QMessageBox::warning(nullptr, u8"警告", u8"发现 configuration Name:" + it.key() + u8" , 发现问题，请检查"); + return false; + } + } + } + + + if (!this->checkPRFModel() || !this->resizePRFPluse()) { + return false; + } + + return true; + +} + +void LAMPTOOLAPI FEKOBase::NearFieldEchoCSVParser::toThetapolar(const QString& filePath) +{ + + this->toEchoData(filePath, 2); +} + +void LAMPTOOLAPI FEKOBase::NearFieldEchoCSVParser::toPhiPolar(const QString& filePath) +{ + this->toEchoData(filePath, 1); +} + +void LAMPTOOLAPI FEKOBase::NearFieldEchoCSVParser::toRPolar(const QString& filePath) +{ + this->toEchoData(filePath, 0); +} + +void LAMPTOOLAPI FEKOBase::NearFieldEchoCSVParser::saveCSV(const QString& filePath) +{ +} + +/// +/// 输出脉冲回波文件 +/// +/// 文件路径 +/// 0: R,1: phi,2: theta +void LAMPTOOLAPI FEKOBase::NearFieldEchoCSVParser::toEchoData(const QString& filePath, size_t outDataName) +{ + // 脉冲整理 + Eigen::MatrixXcd echoData = Eigen::MatrixXcd::Zero(this->prfData.size(), this->freqPoints); + Eigen::MatrixXd antpos = Eigen::MatrixXd::Zero(this->prfData.size(), 5); //x,y,z,theta,phi" + + for (size_t i = 0; i < this->prfData.size(); i++) { + antpos(i, 0) = this->prfData[i].px; + antpos(i, 1) = this->prfData[i].py; + antpos(i, 2) = this->prfData[i].pz; + antpos(i, 3) = this->prfData[i].electricFieldDataList[0].theta; + antpos(i, 4) = this->prfData[i].electricFieldDataList[0].phi; + for (size_t j = 0; j < this->freqPoints; j++) { + if (outDataName == 0) { + echoData(i, j) = std::complex(this->prfData[i].electricFieldDataList[j].reEr, this->prfData[i].electricFieldDataList[j].imEr); + } + else if (outDataName == 1) { + echoData(i, j) = std::complex(this->prfData[i].electricFieldDataList[j].reEphi, this->prfData[i].electricFieldDataList[j].imEphi); + } + else { + echoData(i, j) = std::complex(this->prfData[i].electricFieldDataList[j].reEtheta, this->prfData[i].electricFieldDataList[j].imEtheta); + } + + } + } + + EchoDataClass echodataTemp; + echodataTemp.setAntPos(antpos); + echodataTemp.setFreqStart(this->freqStart); + echodataTemp.setFreqEnd(this->freqEnd); + echodataTemp.setFreqpoints(this->freqPoints); + echodataTemp.setEchoData(echoData); + echodataTemp.SaveEchoData(filePath); +} + +LAMPTOOLAPI FEKOBase::EchoDataClass::EchoDataClass(const FEKOBase::EchoDataClass& inecho) +{ + + this->echoData = inecho.getEchoData(); + this->antPos = inecho.getAntPos(); + this->freqStart = inecho.getFreqStart(); + this->freqEnd = inecho.getFreqEnd(); + this->freqpoints = inecho.getFreqpoints(); +} + + + +LAMPTOOLAPI FEKOBase::EchoDataClass::EchoDataClass() +{ + +} + +LAMPTOOLAPI FEKOBase::EchoDataClass::~EchoDataClass() +{ +} + + + + +void LAMPTOOLAPI FEKOBase::EchoDataClass::setEchoData(Eigen::MatrixXcd echoData) +{ + this->echoData = echoData; +} + +Eigen::MatrixXcd LAMPTOOLAPI FEKOBase::EchoDataClass::getEchoData() const +{ + return this->echoData; +} + +void LAMPTOOLAPI FEKOBase::EchoDataClass::setAntPos(Eigen::MatrixXd antPos) +{ + this->antPos = antPos; +} + +Eigen::MatrixXd LAMPTOOLAPI FEKOBase::EchoDataClass::getAntPos() const +{ + return this->antPos; +} + +void LAMPTOOLAPI FEKOBase::EchoDataClass::setFreqStart(double freqStart) +{ + this->freqStart = freqStart; +} + +double LAMPTOOLAPI FEKOBase::EchoDataClass::getFreqStart() const +{ + return this->freqStart; +} + +void LAMPTOOLAPI FEKOBase::EchoDataClass::setFreqEnd(double freqEnd) +{ + this->freqEnd = freqEnd; +} + +double LAMPTOOLAPI FEKOBase::EchoDataClass::getFreqEnd() const +{ + return this->freqEnd; +} + +void LAMPTOOLAPI FEKOBase::EchoDataClass::setFreqpoints(int freqpoints) +{ + this->freqpoints = freqpoints; +} + +int LAMPTOOLAPI FEKOBase::EchoDataClass::getFreqpoints() const +{ + return this->freqpoints; +} + +void LAMPTOOLAPI FEKOBase::EchoDataClass::loadEchoData(const QString& filePath) +{ + std::ifstream file(reinterpret_cast(filePath.utf16()), std::ios::binary); + + if (file.is_open()) { + // 保存频率变量 + INT32 PRFRow = 0; + INT32 freqCol = 0; + // 读取其他成员变量 + file.read(reinterpret_cast(&freqStart), sizeof(double)); + file.read(reinterpret_cast(&freqEnd), sizeof(double)); + file.read(reinterpret_cast(&freqpoints), sizeof(INT32)); + file.read(reinterpret_cast(&PRFRow), sizeof(INT32)); + file.read(reinterpret_cast(&freqCol), sizeof(INT32)); + echoData = Eigen::MatrixXcd::Zero(PRFRow, freqCol); // 重新分配内存 + antPos = Eigen::MatrixXd::Zero(PRFRow, 5); + file.read(reinterpret_cast(antPos.data()), sizeof(double) * PRFRow * 5); // 保存天线位置 + file.read(reinterpret_cast(echoData.data()), sizeof(std::complex) * PRFRow * freqCol); // 保存回波数据 + file.close(); // 关闭文件 + + this->freqStart = freqStart; + this->freqEnd = freqEnd; + this->freqpoints = freqpoints; + this->echoData = echoData; + this->antPos = antPos; + qDebug() << u8"回波数据加载完成"; + } + else { + qWarning() << "Error: Unable to open file for reading."; + } +} + +void LAMPTOOLAPI FEKOBase::EchoDataClass::SaveEchoData(const QString& filePath) +{ + if (echoData.rows() != antPos.rows() || antPos.cols() != 5) { + qDebug() << "Error: antPos.size()!=echoData.rows()*5" + + QString::number(echoData.rows()) + u8"!=" + QString::number(antPos.rows()) + + QString::number(antPos.cols()) + u8"!=5"; + return; + } + else {} + std::ofstream file(reinterpret_cast(filePath.utf16()), std::ios::binary); + qDebug() << "==================================================================="; + qDebug() << u8"回波脉冲文件格式说明："; + qDebug() << u8"|freqStart(double)，freqEnd(double)，freqpoints(int)|PRFROW(int),PRFCol(int)|antPos|echodata|"; + qDebug() << u8"其中antPos: PRFrow*3 double x,y,z,incidence,azangle,theta,phi"; + qDebug() << u8"其中echodata: PRFrow*PRFCol complex "; + qDebug() << "==================================================================="; + if (file.is_open()) { + INT32 PRFRow = echoData.rows(); + INT32 freqCol = echoData.cols(); + file.write(reinterpret_cast(&freqStart), sizeof(double)); // 保存频率变量 + file.write(reinterpret_cast(&freqEnd), sizeof(double)); + file.write(reinterpret_cast(&freqpoints), sizeof(INT32)); + file.write(reinterpret_cast(&PRFRow), sizeof(INT32)); + file.write(reinterpret_cast(&freqCol), sizeof(INT32)); + file.write(reinterpret_cast(antPos.data()), sizeof(double) * antPos.size()); // 保存天线位置 + file.write(reinterpret_cast(echoData.data()), sizeof(std::complex) * echoData.size()); // 保存回波数据 + file.close(); // 关闭文件 + + } + else { + qDebug() << "Error: Unable to open file for writing."; + } +} + + + +Eigen::MatrixXd LAMPTOOLAPI FEKOBase::WINDOWFun(Eigen::MatrixXcd& echo, ImageAlgWindowFun winfun) +{ + size_t Nxa = echo.rows(); + size_t Nf = echo.cols(); + Eigen::MatrixXd normw = Eigen::MatrixXd::Zero(Nxa, Nf); // 归一化窗函数 + + if (winfun == ImageAlgWindowFun::HANMMING) { + normw = Hanning(Nf, Nxa); + echo = echo.array() * normw.array(); + } + else { + normw = normw.array() * 0 + 1 / (Nxa * Nf); + } + + return normw; +} + +QList LAMPTOOLAPI FEKOBase::getFEKOImageAlgorithmList() +{ + QList list; + FEKOImageAlgorithm alg; + for (alg = FEKOImageAlgorithm::TBP_TIME; alg < FEKOImageAlgorithm::UNKONW; alg = (FEKOImageAlgorithm)(alg + 1)) { + list.append(FEKOImageAlgorithm2String(alg)); + } + return list; +} + +FEKOBase::FEKOImageAlgorithm LAMPTOOLAPI FEKOBase::String2FEKOImageAlgorithm(QString str) +{ + FEKOImageAlgorithm alg; + for (alg = FEKOImageAlgorithm::TBP_TIME; alg < FEKOImageAlgorithm::UNKONW; alg = (FEKOImageAlgorithm)(alg + 1)) { + if (str == FEKOImageAlgorithm2String(alg)) { + return alg; + } + else {} + } + return alg; +} + +QString LAMPTOOLAPI FEKOBase::FEKOImageAlgorithm2String(FEKOImageAlgorithm alg) +{ + // 将FEKOImageAlgorithm 枚举转换为字符串 + switch (alg) + { + case FEKOBase::TBP_TIME: + return u8"TBP_TIME"; + case FEKOBase::TBP_FREQ: + return u8"TBP_FREQ"; + default: + return u8"UNKONW"; + } +} + +// 请仿照FEKOImageAlgorithm枚举的写法，构建QString 与 ImageAlgWindowFun 的转换函数 + +QList LAMPTOOLAPI FEKOBase::getImageAlgWindowFunList() +{ + QList list; + ImageAlgWindowFun alg; + for (alg = ImageAlgWindowFun::NOWINDOWS; alg < ImageAlgWindowFun::UNKONWWINDOW; alg = (ImageAlgWindowFun)(alg + 1)) { + list.append(ImageAlgWindowFun2String(alg)); + } + return list; +} + +FEKOBase::ImageAlgWindowFun LAMPTOOLAPI FEKOBase::String2ImageAlgWindowFun(QString str) +{ + ImageAlgWindowFun alg; + for (alg = ImageAlgWindowFun::UNKONWWINDOW; alg < ImageAlgWindowFun::UNKONWWINDOW; alg = (ImageAlgWindowFun)(alg + 1)) { + if (str == ImageAlgWindowFun2String(alg)) { + return alg; + } + else {} + } + return alg; +} + +QString LAMPTOOLAPI FEKOBase::ImageAlgWindowFun2String(FEKOBase::ImageAlgWindowFun alg) +{ + switch (alg) + { + case FEKOBase::NOWINDOWS: + return u8"NOWINDOWS"; + case FEKOBase::HANMMING: + return u8"HANMMING"; + default: + return u8"UNKONWWINDOW"; + } +} + + +bool LAMPTOOLAPI FEKOBase::BPImage_TIME(QString& restiffpath, Eigen::MatrixXcd& echoData, Eigen::MatrixXd& antPos, Eigen::MatrixXd& freqmatrix, Eigen::MatrixXd& X, Eigen::MatrixXd& Y, Eigen::MatrixXd& Z, ImageAlgWindowFun winfun) +{ + // BP成像算法 + const double c = 0.299792458; // 光速 + + return true; +} + +bool LAMPTOOLAPI FEKOBase::FBPImage_FREQ(QString& restiffpath, Eigen::MatrixXcd& echoData, Eigen::MatrixXd& antPos, Eigen::MatrixXd& freqmatrix, Eigen::MatrixXd& X, Eigen::MatrixXd& Y, Eigen::MatrixXd& Z, ImageAlgWindowFun winfun) +{ + + // BP成像算法 + const double c = 0.299792458; // 光速 + const std::complex j(0, 1); // 虚数单位 + const size_t image_height = X.rows(); + const size_t image_width = X.cols(); + const size_t PRFCount = echoData.rows(); + const size_t frepoints = echoData.cols(); + Eigen::MatrixXcd factorj = Eigen::MatrixXcd::Zero(1, echoData.cols()); // 1 , freqs + factorj = freqmatrix.array().cast>().array() * j * 4 * M_PI / c; // 声明校正矩阵 + Eigen::MatrixXcd im_final = Eigen::MatrixXcd::Zero(image_height, image_width); + + + +#ifdef __SHOWPROCESS // 进度条展示 + QProgressDialog progressDialog(u8"BP_FREQ成像进度", u8"终止", 0, image_height); + progressDialog.setWindowTitle(u8"成像中"); + progressDialog.setWindowModality(Qt::WindowModal); + progressDialog.setAutoClose(true); + progressDialog.setValue(0); + progressDialog.setMaximum(image_height); + progressDialog.setMinimum(0); + progressDialog.show(); +#endif // __SHOWPROCESS + + + +#ifdef __SHOWIMAGEPROCESSRESULT // 加窗处理 + Eigen::MatrixXd normw = WINDOWFun(echoData, winfun); +#endif + +#ifdef __IMAGEPARALLEL +#pragma omp parallel for + for (long ii = 0; ii < image_height; ii++) + { + Eigen::MatrixXd im_R = Eigen::MatrixXd::Zero(PRFCount, 1);// 图像到脉冲矩阵 prf*1 + Eigen::MatrixXcd term_R = Eigen::MatrixXcd::Zero(1, frepoints); + Eigen::MatrixXcd pluse_R = Eigen::MatrixXcd::Zero(PRFCount, 1); + + for (size_t jj = 0; jj < image_width; jj++) { + im_R.col(0) = ((antPos.col(0).array() - X(ii, jj)).array().pow(2) + + (antPos.col(1).array() - Y(ii, jj)).array().pow(2) + + (antPos.col(2).array() - Z(ii, jj)).array().pow(2)) + .array().sqrt().array(); // 获取目标到天线矩阵 计算 ok + // im_R PRF_count x 1 + for (size_t tt = 0; tt < PRFCount; tt++) { + term_R = (im_R(tt, 0) * factorj.array()).array().exp(); // ok + pluse_R(tt, 0) = (echoData.row(tt).array() * term_R.array()).sum(); // 对每个脉冲校正求和 + } + im_final(ii, jj) = pluse_R.array().sum(); // 计算图像 + } + } +#else + + Eigen::MatrixXd im_R = Eigen::MatrixXd::Zero(PRFCount, 1);// 图像到脉冲矩阵 prf*1 + Eigen::MatrixXcd term_R = Eigen::MatrixXcd::Zero(1, frepoints); + Eigen::MatrixXcd pluse_R = Eigen::MatrixXcd::Zero(PRFCount, 1); + for (long ii = 0; ii < image_height; ii++) + { + for (size_t jj = 0; jj < image_width; jj++) { + im_R.col(0) = ((antPos.col(0).array() - X(ii, jj)).array().pow(2) + + (antPos.col(1).array() - Y(ii, jj)).array().pow(2) + + (antPos.col(2).array() - Z(ii, jj)).array().pow(2)) + .array().sqrt().array(); // 获取目标到天线矩阵 计算 ok + // im_R PRF_count x 1 + for (size_t tt = 0; tt < PRFCount; tt++) { + term_R = (im_R(tt, 0) * factorj.array()).array().exp(); // ok + pluse_R(tt, 0) = (echoData.row(tt).array() * term_R.array()).sum(); // 对每个脉冲校正求和 + } + im_final(ii, jj) = pluse_R.array().sum(); // 计算图像 + } +#ifdef __SHOWPROCESS + progressDialog.setValue(ii); +#endif + } +#endif + +#ifdef __SHOWIMAGEPROCESSRESULT // 加窗处理 + im_final = im_final.array() / (normw.array().sum()); + +#endif + + + // 3. 保存图像 + + + QString newrestiffpath = restiffpath; + + + Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2,3); + gt(0, 0) = X(0, 0); + gt(0, 1) = X(0, 1) - X(0, 0); + gt(0, 2) = 0; + gt(1, 0) = Y(0, 0); + gt(1, 1) = 0; + gt(1, 2) = Y(0, 0) - Y(1, 0); + std::cout << gt << std::endl; + gdalImageComplex im_final_gdal = CreategdalImageComplex(restiffpath,im_final.rows(),im_final.cols(),1,gt, QString(u8""), true,true); + im_final_gdal.setData(im_final); + im_final_gdal.saveImage(); + + qDebug() << QString(u8"保存文件地址：" + restiffpath); +#ifdef __SHOWPROCESS + progressDialog.setWindowTitle(u8"保存文件：" + restiffpath); + progressDialog.setValue(image_height); + progressDialog.close(); +#endif + + // 绘制图像 +#ifdef __SHOWIMAGEPROCESSRESULT + // 将复数矩阵的幅度转换为实数矩阵 + matplot::figure(); + matplot::vector_2d amplitude_im_final(im_final.rows(), matplot::vector_1d(im_final.cols())); + for (int i = 0; i < im_final.rows(); ++i) { + for (int j = 0; j < im_final.cols(); ++j) { + amplitude_im_final[i][j] = std::abs(im_final(i, j)); + } + } + matplot::imagesc(amplitude_im_final); + matplot::colorbar(); + matplot::show(); +#endif // __SHOWMATPLOTLIBCPP + + return true; +} + + + diff --git a/src/LAMPTool/FEKOBaseToolClass.h b/src/LAMPTool/FEKOBaseToolClass.h new file mode 100644 index 0000000..ddbe460 --- /dev/null +++ b/src/LAMPTool/FEKOBaseToolClass.h @@ -0,0 +1,344 @@ +#pragma once +#ifndef _FEKOBASETOOLCLASS_H_ +#define _FEKOBASETOOLCLASS_H_ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "referenceHeader.h" +#include "GeoOperator.h" + + +/** +* 此类用于FEKO的基本类型，主要功能如下： +* 1. 根据分辨率、中心频率、带宽，近远斜距，计算频率起始点，以及采样点数 +* 2. 根据起始成像点、成像时长、PRF采样点数，获取PRF采样点数 +* 3. 结合 PRF采样点坐标，入射角，以及与飞行方向的夹角（默认为 90度）, 参考高程，计算入射波的参考平面部分 +* 4. +***/ + +namespace FEKOBase { + + //========================================================== + // FEKO常用坐标系 + //========================================================== + enum FEKOCoordinateSystem { + Spherical, // 球坐标系 + Cartesian, // 笛卡尔坐标系 + UNKONWFEKOCOORDINATESYSTEM // 必须为最后一个表示未知 + }; + FEKOBase::FEKOCoordinateSystem LAMPTOOLAPI FEKOCoordinateSystemString2Enum(QString str); + QString LAMPTOOLAPI QString2FEKOCoordinateSystem(FEKOBase::FEKOCoordinateSystem mode); + + + //========================================================== + // FEKO成像模式枚举 + //========================================================== + enum FEKOImageMode + { + Strip, + Scane, + ISAR, + CircleSAR, + UNKNOW + }; + + FEKOImageMode LAMPTOOLAPI FEKOImageModeString2Enum(QString str); + QString LAMPTOOLAPI FEKOImageModeenumToString(FEKOImageMode mode); + + + + + //========================================================== + // FEKO成像仿真参数类，主要用来搭建统一的仿真成像模块 + //========================================================== + + + + + //========================================================== + // 频率参数 + //========================================================== + struct freqParams { // 频率参数 + double startfreqs; + double endfreqs; + size_t freqpoint; + }; + + //========================================================== + // 输入的卫星参数 + //========================================================== + struct SatellitePosition { // 卫星姿态 + double Px = 0, Py = 0, Pz = 0; + }; + struct SatelliteVelocity { // 卫星速度 + double Vx = 0, Vy = 0, Vz = 0; + }; + + struct SatelliteState { // 卫星矢量 + SatellitePosition pos; + SatelliteVelocity vel; + }; + + + //========================================================== + // FEKO 远场等效源姿态参数 + //========================================================== + struct FEKOantPitionDirect { + double x = 0; + double y = 0; + double z = 0; + double theta = 0; + double phi = 0; + }; + + + //=========================================================== + // FEKO成像设置参数 + //=========================================================== + struct FEKOImageSettingParams { + double min_x = 0; + double max_x = 0; + double min_y = 0; + double max_y = 0; + double plane_z = 0; // 平面高程 + size_t ImageWidth = 0; + size_t ImageHeight = 0; + + }; + + + + //========================================================== + // FEKO参数文件等效参数 + //========================================================== + struct FEKOSatelliteParams { // FEKOPRF脉冲参数 + size_t PRFidx = 0; // PRF 脉冲计数 + SatelliteState pose;// 卫星矢量 + double incidenceAngle = 0; // 入射角 + double AzAngle = 0; // 姿态角 + FEKOantPitionDirect antpos;// 天线实际姿态 + bool isRight = false; // 判断左右视，false 左，right 右 + }; + + + ///////////////////////////////////////////////// + ///// 函数类 + ////////////////////////////////////////////// + /// + /// 频率Setting + /// + /// 中心频率 GHZ + /// 分辨率 米 + /// 带宽 GHz + /// 分辨率 米 + /// + /// + freqParams LAMPTOOLAPI getFreqSetting(double centerFreq, double resolution, double bandWidth, double scenceRange, bool isResolution = false); + + FEKOSatelliteParams LAMPTOOLAPI createFEKOSatelliteParams(double Px, double Py, double Pz, double Vx, double Vy, double Vz, double incidenceAngle, double AzAngle, double theta, double phi, bool isRight, size_t PRFIdx = 0); + FEKOSatelliteParams LAMPTOOLAPI createFEKOSatelliteParams(SatelliteState pose, double incidenceAngle, double AzAngle, FEKOantPitionDirect antpos, size_t PRFIdx = 0); + + SatelliteState LAMPTOOLAPI FEKOSatelliteParams2SatelliteState(FEKOSatelliteParams parmas); + FEKOantPitionDirect LAMPTOOLAPI FEKOSatelliteParams2FEKOantPitionDirect(FEKOSatelliteParams parmas); + + /// + /// 将卫星姿态转换为 FEKO 天线坐标， + /// 注意默认模型的Z轴为天线的指向 + /// + /// 卫星矢量 + /// 入射角 + /// 侧视角 + /// 是否为右视 + /// 结果文件 + /// 输入模型 + /// + TopoDS_Shape LAMPTOOLAPI SatellitePos2FEKOAntPos(SatelliteState satepos, double incidenceAngle, double AzAngle, bool isRIGHT, FEKOantPitionDirect* antposition_Direct, TopoDS_Shape inDs); + + + + //=============================================== + // FEKO结果解析文件 + //=============================================== + + struct ElectricFieldData { + QString fileType; + QString fileFormat; + QString source; + QString date; + double radius; + double theta; + double phi; + double reEr; + double imEr; + double reEtheta; + double imEtheta; + double reEphi; + double imEphi; + QString configurationName; + QString requestName; + double frequency; + QString coordinateSystem; + Point_3d origin; + int numRadiusSamples; + int numThetaSamples; + int numPhiSamples; + QString resultType; + int numHeaderLines; + size_t prfidx = -1; // 脉冲计数，>0 + }; + + struct PRFPluseData { // 单个PRF脉冲数据格式 + size_t prfidx; //脉冲次数 + double freqstart; + double freqend; + size_t freqpoints; + + double px, py, pz; + //double theta, phi; + //double incidence, azangle; + //std::vector freqlist; + std::vector electricFieldDataList; // 单频点信息 + }; + + bool LAMPTOOLAPI compareElectricFieldDataInFreq(const ElectricFieldData& a, const ElectricFieldData& b); + bool LAMPTOOLAPI comparePRFPluseDataInPRFIdx(const PRFPluseData& a, const PRFPluseData& b); + + class LAMPTOOLAPI NearFieldEchoCSVParser { + public: + NearFieldEchoCSVParser(); + ~NearFieldEchoCSVParser(); + private: + bool usePRFCountMode = true; + //std::vector electricFieldDataList; + QMap> electricFieldDataList; // 电场数据 + //QMap prfPluseMap; + std::vector prfData; + size_t freqPoints; + double freqStart; + double freqEnd; + // 频率参数 + private: // 内部检查函数 + bool checkPRFModel(); + bool resizePRFPluse(); //回波整理 + + public: + bool parseCSV(const QString& filePath); // 读取回波数据文件 + void toThetapolar(const QString& filePath);// 输出theta 极化 + void toPhiPolar(const QString& filePath);// 输出phi 极化 + void toRPolar(const QString& filePath);// 输出phi 极化 + void saveCSV(const QString& filePath);// 输出csv文件 + private: + void toEchoData(const QString& filePath, size_t outDataName);// 输出回波数据 + }; + + + + + //======================================================================== + // 成像回波格式 + // file type: + // freqStart,freqEnd,freqPoint,isRight, + // PRF1,Pos,incidenceAngle,AzAngle,echoDatalist + // PRF2,Pos,incidenceAngle,AzAngle,echoDatalist + // 。 + // 。 + // 。 + // 注意Bp并不关心脉冲的顺序，只是关注脉冲的坐标位置，默认按照脉冲的解析顺序进行组织 + //======================================================================== + class LAMPTOOLAPI EchoDataClass { + private: // 成像变量 + Eigen::MatrixXcd echoData; // 回波数据 + Eigen::MatrixXd antPos;// 每个脉冲的坐标 + double freqStart; // 起始频率 + double freqEnd; // 终止频率 + int freqpoints; // 频率点数 + + public: + EchoDataClass(const FEKOBase::EchoDataClass& inecho); + EchoDataClass(); + ~EchoDataClass(); + public: + // 根据每个成员变量构建属性 + void setEchoData(Eigen::MatrixXcd echoData); + Eigen::MatrixXcd getEchoData() const; + void setAntPos(Eigen::MatrixXd antPos); + Eigen::MatrixXd getAntPos() const; + void setFreqStart(double freqStart); + double getFreqStart() const; + void setFreqEnd(double freqEnd); + double getFreqEnd() const; + void setFreqpoints(int freqpoints); + int getFreqpoints() const; + + void loadEchoData(const QString& filePath); // 加载回波数据文件 + void SaveEchoData(const QString& filePath); // 保存回波数据文件 + }; + + + + + + //========================================================== + // 仿真成像算法类 + // BP成像算法 + // 建议将所有的频率处理到 GHz 单位，避免频率过大导致的精度问题 + // 回波矩阵： + // PRF1:f1,f2,f3,f4 + // PRF2：f1,f2,f3,f4 + // PRF2：f1,f2,f3,f4 + //========================================================== + + enum ImageAlgWindowFun // 成像方法加窗方法 + { + NOWINDOWS, + HANMMING, + UNKONWWINDOW + }; + + /// + /// 加窗 + /// + /// 行:脉冲，列：频点 + /// + Eigen::MatrixXd LAMPTOOLAPI WINDOWFun(Eigen::MatrixXcd& echo, ImageAlgWindowFun winfun = ImageAlgWindowFun::HANMMING); + + enum FEKOImageAlgorithm + { + TBP_TIME, + TBP_FREQ, + UNKONW // 必须为最后一个表示未知 + }; + + QList LAMPTOOLAPI getFEKOImageAlgorithmList(); + FEKOImageAlgorithm LAMPTOOLAPI String2FEKOImageAlgorithm(QString str); + QString LAMPTOOLAPI FEKOImageAlgorithm2String(FEKOImageAlgorithm alg); + + // 请仿照FEKOImageAlgorithm枚举的写法，构建QString 与 ImageAlgWindowFun 的转换函数 + QList LAMPTOOLAPI getImageAlgWindowFunList(); + ImageAlgWindowFun LAMPTOOLAPI String2ImageAlgWindowFun(QString str); + QString LAMPTOOLAPI ImageAlgWindowFun2String(ImageAlgWindowFun alg); + + + bool LAMPTOOLAPI BPImage_TIME(QString& restiffpath, Eigen::MatrixXcd& echoData, Eigen::MatrixXd& antPos, Eigen::MatrixXd& freqmatrix, Eigen::MatrixXd& X, Eigen::MatrixXd& Y, Eigen::MatrixXd& Z, ImageAlgWindowFun winfun = ImageAlgWindowFun::HANMMING); // BP成像 + bool LAMPTOOLAPI FBPImage_FREQ(QString& restiffpath, Eigen::MatrixXcd& echoData, Eigen::MatrixXd& antPos, Eigen::MatrixXd& freqmatrix, Eigen::MatrixXd& X, Eigen::MatrixXd& Y, Eigen::MatrixXd& Z, ImageAlgWindowFun winfun = ImageAlgWindowFun::HANMMING); // FBP成像 + + + +} + + +#endif \ No newline at end of file diff --git a/src/LAMPTool/FEKONearBPBasic.cpp b/src/LAMPTool/FEKONearBPBasic.cpp new file mode 100644 index 0000000..d373bb6 --- /dev/null +++ b/src/LAMPTool/FEKONearBPBasic.cpp @@ -0,0 +1,1028 @@ +#include "FEKONearBPBasic.h" +#include "interpolation.h" +#include "FileOperator.h" +#include "BaseTool.h" +#include "ImageOperatorBase.h" +#include "SARImageBase.h" +#include +#include +#include +#include +#include +#include + +#include "SARBaseTool.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + + +// 定义插值函数，以处理复数值 +std::complex LAMPTOOLAPI InterpolateComplex(Eigen::MatrixXd& xi, Eigen::MatrixXd& yi, Eigen::MatrixXcd& data, double x, double y) { + int m = xi.rows(); + int n = xi.cols(); + + if (x < xi(0, 0) || x > xi(m - 1, n - 1) || y < yi(0, 0) || y > yi(m - 1, n - 1)) { + std::complex defaultMatrix = std::complex(0,0); // 超出范围时返回默认值 + return defaultMatrix; + } + + int i, j; + for (i = 0; i < m - 1; i++) { + if (x >= xi(i, 0) && x <= xi(i + 1, 0)) { + break; + } + } + + for (j = 0; j < n - 1; j++) { + if (y >= yi(0, j) && y <= yi(0, j + 1)) { + break; + } + } + + double x1 = xi(i, 0); + double x2 = xi(i + 1, 0); + double y1 = yi(0, j); + double y2 = yi(0, j + 1); + + std::complex Q11 = data(i, j); + std::complex Q12 = data(i, j + 1); + std::complex Q21 = data(i + 1, j); + std::complex Q22 = data(i + 1, j + 1); + + double dx1 = x2 - x; + double dx2 = x - x1; + double dy1 = y2 - y; + double dy2 = y - y1; + double denom = (x2 - x1) * (y2 - y1); + + std::complex interpolatedValue = (Q11 * dx1 * dy1 + Q21 * dx2 * dy1 + Q12 * dx1 * dy2 + Q22 * dx2 * dy2) / denom; + return interpolatedValue; +} + + +Eigen::MatrixXd LAMPTOOLAPI Cartesian2Spherical(Eigen::MatrixXd CartesianPoint) +{ + Eigen::MatrixXd result = CartesianPoint; + size_t rows = result.rows(); + size_t cols = result.cols(); + for (int i = 0; i < rows; i++) { + double x = CartesianPoint(i, 0); + double y = CartesianPoint(i, 1); + double z = CartesianPoint(i, 2); + result(i, 0) = std::sqrt(CartesianPoint(i, 0) * CartesianPoint(i, 0) + + CartesianPoint(i, 1) * CartesianPoint(i, 1) + + CartesianPoint(i, 2) * CartesianPoint(i, 2)); + result(i, 1) = std::acos(CartesianPoint(i, 2) / result(i, 0));// 行 theta + if (CartesianPoint(i, 0) == 0 && CartesianPoint(i, 1) == 0) { + result(i, 2) = 0; + } + else { + result(i, 2) = std::asin(y / (std::sin(result(i, 1)) * result(i, 0)));//一 三 -pi/2 pi/2 + if (CartesianPoint(i, 0) < 0 && CartesianPoint(i, 1) > 0) { // 二 + result(i, 2) = PI - result(i, 2);// phi + } + else if (CartesianPoint(i, 0) < 0 && CartesianPoint(i, 1) < 0) { + result(i, 2) = -1 * PI - result(i, 2);// phi + } + else {} + } + } + + return result; +} + +int LAMPTOOLAPI BP2DProcess(QString in_path, QString out_path, double Rref, double minX, double maxX, double minY, double maxY, double PlaneZ, int ImageHeight, int ImageWidth) +{ + BP2DProcessClass process; + process.initProcess(in_path, out_path, Rref, minX, maxX, minY, maxY, PlaneZ, ImageHeight, ImageWidth); + process.start(); + return -1; +} + + +int LAMPTOOLAPI FBP2DProcess(QString in_path, QString out_path, double Rref, double minX, double maxX, double minY, double maxY, double PlaneZ, int ImageHeight, int ImageWidth) +{ + FBP2DProcessClass process; + process.initProcess(in_path, out_path, Rref, minX, maxX, minY, maxY, PlaneZ, ImageHeight, ImageWidth); + process.start(); + return -1; +} + +int LAMPTOOLAPI build2Bin(QString path, int width, int height, std::vector& freqs, Eigen::MatrixXd AntPostion, Eigen::MatrixXcd echo) +{ + /* +* % 构建回波结果 +fid=fopen("H_echo.bin","w+"); +fwrite(fid,801,'int32'); +fwrite(fid,801,'int32'); +for i=1:length(freqs) + fwrite(fid,freqs(i),'double'); +end + +for i=1:length(AntPostion) + fwrite(fid,AntPostion(i,1),'double'); + fwrite(fid,AntPostion(i,2),'double'); + fwrite(fid,AntPostion(i,3),'double'); + fwrite(fid,real(echo_H_datas(i,:)),'double'); + fwrite(fid,imag(echo_H_datas(i,:)),'double'); +end +fclose(fid); +**/ +// 构建二进制文件 + + std::ofstream file(path.toUtf8().constData(), std::ios::binary); + double freq; + + file.write(reinterpret_cast(&height), sizeof(height)); + file.write(reinterpret_cast(&width), sizeof(width)); + + for (size_t i = 0; i < freqs.size(); i++) { + freq = freqs[i]; + file.write(reinterpret_cast(&freq), sizeof(freq)); + } + + double x, y, z; + double real_val, imag_val; + + for (size_t i = 0; i < AntPostion.rows(); i++) { + x = AntPostion(i, 0); + y = AntPostion(i, 1); + z = AntPostion(i, 2); + file.write(reinterpret_cast(&x), sizeof(x)); + file.write(reinterpret_cast(&y), sizeof(y)); + file.write(reinterpret_cast(&z), sizeof(z)); + for (size_t j = 0; j < echo.cols(); j++) { + real_val = std::real(echo(i, j)); + imag_val = std::imag(echo(i, j)); + file.write(reinterpret_cast(&real_val), sizeof(real_val)); + file.write(reinterpret_cast(&imag_val), sizeof(imag_val)); + } + } + + file.close(); + return 0; +} + +template +Eigen::MatrixXcd LAMPTOOLAPI BP2DImageByPluse(Eigen::MatrixXcd timeEcho, Eigen::MatrixXd AntPosition, double minX, double maxX, double minY, double maxY, double PlaneZ, double Rref, size_t ImageWidth, size_t ImageHeight, double startfreq, size_t timefreqnum, double timeFreqBandWidth, T* logclss) +{ + bool logfun = !(nullptr == logclss); // 空指针 + double delta_x = (maxX - minX) / (ImageWidth - 1); + double delta_y = (maxY - minY) / (ImageHeight - 1); + + double* Rs = new double[timefreqnum]; // 时域 距离插值 + for (int i = 0; i < timefreqnum; i++) { + Rs[i] = (LIGHTSPEED / timeFreqBandWidth) * i; // 参考 《雷达成像算法》-- 匹配滤波一章 + } + double wave_len = LIGHTSPEED / startfreq; + if (logfun) { + logclss->logFUN(0, "process....."); + } + + size_t process = 0; + size_t PRFNUM = timeEcho.rows(); + Eigen::MatrixXcd img = Eigen::MatrixXcd::Zero(ImageHeight, ImageWidth); + + for (int i = 0; i < PRFNUM; i++) { + gsl_interp_accel* accopm = gsl_interp_accel_alloc(); + const gsl_interp_type* tt = gsl_interp_linear; //gsl_interp_cspline_periodic; + gsl_spline* spline_time_real = gsl_spline_alloc(tt, timefreqnum); + gsl_spline* spline_time_imag = gsl_spline_alloc(tt, timefreqnum); + + double* real_time_ys = new double[timefreqnum]; + double* imag_time_ys = new double[timefreqnum]; + for (int tt = 0; tt < timefreqnum; tt++) { + real_time_ys[tt] = timeEcho(i, tt).real(); + imag_time_ys[tt] = timeEcho(i, tt).imag(); + } + gsl_spline_init(spline_time_real, Rs, real_time_ys, timefreqnum); + gsl_spline_init(spline_time_imag, Rs, imag_time_ys, timefreqnum); + for (int ii = 0; ii < ImageHeight; ii++) { + for (int jj = 0; jj < ImageWidth; jj++) { + double px = ii * delta_x + minX; + double py = jj * delta_y + minY; + double pz = PlaneZ; + double dR = 0; + double R = std::sqrt( + std::pow(AntPosition(i, 0) - px, 2) + + std::pow(AntPosition(i, 1) - py, 2) + + std::pow(AntPosition(i, 2) - pz, 2) + ); + dR = R - Rref; + std::complex echo_temp( + gsl_spline_eval(spline_time_real, R, accopm),// 插值实部 + gsl_spline_eval(spline_time_imag, R, accopm)// 插值虚部 + ); + img(ii, jj) = img(ii, jj) + echo_temp * std::exp(std::complex(0, 2 * pi * dR / wave_len)); + } + } + gsl_spline_free(spline_time_real); + gsl_spline_free(spline_time_imag); + gsl_interp_accel_free(accopm); + // 保存文件 + process = process + 1; + if (logfun) { + logclss->logFUN((size_t)(process / PRFNUM * 100), "BP process over"); + } + } + if (logfun) { + logclss->logFUN(100, "BP process over"); + } + return img; +} + + +template +Eigen::MatrixXcd LAMPTOOLAPI BP2DImageByPixel(Eigen::MatrixXcd timeEcho, Eigen::MatrixXd AntPosition, double minX, double maxX, double minY, double maxY, double PlaneZ, double Rref, size_t ImageWidth, size_t ImageHeight, double startfreq, size_t timefreqnum, double timeFreqBandWidth, T* logclss) +{ + bool logfun = !(nullptr == logclss); + double delta_x = (maxX - minX) / (ImageWidth - 1); + double delta_y = (maxY - minY) / (ImageHeight - 1); + double maxRange = delta_x * ImageWidth > delta_y * ImageHeight ? delta_x * ImageWidth : delta_y * ImageHeight; + double halfR = Rref * 0.1; + double maxRLs = std::sqrt(Rref * Rref + halfR * halfR) - Rref; + double minR = Rref - maxRange * 2 - maxRLs; // 最小斜距 + double maxR = Rref + maxRange * 2 + maxRLs; // 最大斜距 + + double dftime = timeFreqBandWidth / (timefreqnum - 1); + + // 后续需要增加距离窗的计算方法 + + + double dRs = (LIGHTSPEED / 2.0/timeFreqBandWidth); // 距离向分辨率 + double* Rs = new double[timefreqnum]; // 时域 距离插值 + for (int i = 0; i < timefreqnum; i++) { + Rs[i] = 2*dRs * i; // 参考 《雷达成像算法》-- 匹配滤波一章 + } + //double wave_len = LIGHTSPEED / startfreq; + if (logfun) { + logclss->logFUN(0, "process....."); + } + + + size_t process = 0; + size_t PRFNUM = timeEcho.rows(); + + Eigen::MatrixXcd img = Eigen::MatrixXcd::Zero(ImageHeight, ImageWidth); + omp_lock_t lock; + omp_init_lock(&lock); // 初始化互斥锁 + size_t parallel_step = 100; + + +#pragma omp parallel for // NEW ADD + for (int ii = 0; ii < ImageHeight; ii = ii + parallel_step) { + size_t start_i = ii; + size_t max_i = ii + parallel_step; + max_i = max_i < ImageHeight ? max_i : ImageHeight; + + + double px; + double py ; + double pz; + std::complex echo_sum(0, 0); + double R = 0; + double dR = 0; + + + // 估算坐标的夹角，估算双线性距离 + Eigen::MatrixXd Rline(PRFNUM, 2); + Point_3d Vpluse_p; + Point_3d Vs_e; + double Rline_ref ; + double offer_Rline; + + Eigen::MatrixX> echoline; + double real, imag; + std::complex echo_temp; + std::complex echo_phase(0, 0); + std::complex last_value; + std::complex next_value; + size_t last_ids; + size_t next_ids; + double index = 0; + + Vector3D p; + Vector3D lineP; + Vector3D lineDirection; + double DopplerFreq; + double wave_len; + + + for (int i = start_i; i < max_i; i++) + { + for (int j = 0; j < ImageWidth; j++) { + px = i * delta_x + minX; + py = j * delta_y + minY; + pz = PlaneZ; + echo_sum.real(0); + echo_sum.imag(0); + R = 0; + dR = 0; + + + for (int mm = 0; mm < PRFNUM; mm++) { // 计算每个脉冲情况 + Rline(mm, 0) = std::sqrt( + std::pow(AntPosition(mm, 0) - px, 2) + + std::pow(AntPosition(mm, 1) - py, 2) + + std::pow(AntPosition(mm, 2) - pz, 2) + ); + + // 计算夹角 + Vpluse_p.x = px - AntPosition(mm, 0); + Vpluse_p.y = py - AntPosition(mm, 1); + Vpluse_p.z = pz - AntPosition(mm, 2); + if (mm == 0) { + Vs_e.x = AntPosition(mm + 1, 0) - AntPosition(mm, 0); + Vs_e.y = AntPosition(mm + 1, 1) - AntPosition(mm, 1); + Vs_e.z = AntPosition(mm + 1, 2) - AntPosition(mm, 2); + } + else { + Vs_e.x = AntPosition(mm, 0) - AntPosition(mm - 1, 0); + Vs_e.y = AntPosition(mm, 1) - AntPosition(mm - 1, 1); + Vs_e.z = AntPosition(mm, 2) - AntPosition(mm - 1, 2); + } + + Rline(mm, 1) = Vpluse_p.x * Vs_e.x + Vpluse_p.y * Vs_e.y + Vpluse_p.z * Vs_e.z; // cos(90)=0 + } + Rline_ref = Rref; + // 查找最小值，获取 0多普勒条件下的 最短参考路径 + for (int mm = 0; mm < PRFNUM - 1; mm++) { + if (Rline(mm, 0) == 0) { // 求解最小值 + Rline_ref = Rline(mm, 0); + break; + } + else { + double flag = Rline(mm, 1) * Rline(mm + 1, 1); + if (flag < 0) { // 不同号 + // 求解目标点到 当前直线的最短距离 + p.x = px; + p.y = py; + p.z = pz; + lineP.x = AntPosition(mm, 0); + lineP.y = AntPosition(mm, 1); + lineP.z = AntPosition(mm, 2); + lineDirection.x = AntPosition(mm + 1, 0) - AntPosition(mm, 0); + lineDirection.y = AntPosition(mm + 1, 1) - AntPosition(mm, 1); + lineDirection.z = AntPosition(mm + 1, 2) - AntPosition(mm, 2); + + Rline_ref = pointToLineDistance(p, lineP, lineDirection); + break; + } + } + } + + offer_Rline = Rline_ref - Rref; // 计算偏移线 + Rline = Rline.array() + offer_Rline; + + for (int mm = 0; mm < PRFNUM; mm++) { + R= Rline(mm, 0); // 已经校正了距离徙动线 + + //if (R maxR) { // 距离窗 + // continue; + //} + dR = R - Rref; + // 插值计算 + index = R / 2 / dRs; + last_ids = size_t(std::floor(index)); + next_ids = size_t(std::ceil(index)); + + last_value = timeEcho(mm, last_ids); + next_value = timeEcho(mm, next_ids); + + // 实部，虚部同时插值 + real = last_value.real() + ((next_value.real() - last_value.real()) / (next_ids - last_ids)) * (index - last_ids); + imag = last_value.imag() + ((next_value.imag() - last_value.imag()) / (next_ids - last_ids)) * (index - last_ids); + echo_temp.real(real); + echo_temp.imag(imag); + + DopplerFreq =std::sqrt(R*R - Rline_ref* Rline_ref) / Rline_ref * startfreq; + wave_len = std::abs(LIGHTSPEED / DopplerFreq); + echo_phase.imag(-2 * pi * dR / wave_len); + echo_sum = echo_sum + echo_temp;// *std::exp(echo_phase); + } + img(i, j) = echo_sum; + } + // 保存文件 + + } + omp_set_lock(&lock); + process = process + parallel_step; + process = process < ImageHeight ? process : ImageHeight; + if (logfun) { + logclss->logFUN(size_t(std::round(process * 1.0 / ImageHeight * 100)), "BP process over"); + } + omp_unset_lock(&lock); + } + omp_destroy_lock(&lock); //销毁互斥器 + if (logfun) { + logclss->logFUN(100, "BP process over"); + } + return img; +} + + +int LAMPTOOLAPI BP2DProcessClass::initProcess(QString in_path, QString out_path, double Rref, double minX, double maxX, double minY, double maxY, double PlaneZ, int ImageHeight, int ImageWidth) +{ + // 初始化参数 + this->readEchoFile(in_path); + this->out_path = out_path; + this->Rref = Rref; // 成像中心的参考距离 + this->minX = minX; + this->maxX = maxX; + this->minY = minY; + this->maxY = maxY; + this->Zplane = PlaneZ; + this->ImageHeight = ImageHeight; + this->ImageWidth = ImageWidth; + // 计算坐标的到成像中心的坐标 + this->centerX = (this->maxX + this->minX) / 2; + this->centerY = (this->minY + this->maxY) / 2; + this->AntPosition.col(0) = this->AntPosition.col(0).array() - this->centerX; + this->AntPosition.col(1) = this->AntPosition.col(1).array() - this->centerY; + this->AntPosition.col(2) = this->AntPosition.col(2).array() - this->Zplane; + // 计算频率相关信息 + this->f1 = this->Frequencylist(0); + this->freqnum = this->Frequencylist.rows(); + this->fc = (this->f1 + this->Frequencylist(this->freqnum - 1)); + return -1; +} + +int LAMPTOOLAPI BP2DProcessClass::readEchoFile(QString in_path) +{ + std::ifstream fin(in_path.toUtf8().constData(), std::ios::in | std::ios::binary); + if (!fin.is_open()) exit(2); + // 读取参数 + int height = 0; + int width = 0; + fin.read((char*)&height, sizeof(int)); + fin.read((char*)&width, sizeof(int)); + this->height = height; + this->width = width; + //frequencylist[width*double] + double* freqs = new double[width]; + fin.read((char*)freqs, sizeof(double) * width); + this->Frequencylist = Eigen::VectorXd::Zero(width, 1); // 列向量 + for (int i = 0; i < width; i++) { + this->Frequencylist(i) = freqs[i]; + } + delete[] freqs; + freqs = nullptr; + + // 读取回波并解析脉冲天线位置 + this->AntPosition = Eigen::MatrixXd::Zero(height, 3); + this->echo = Eigen::MatrixXcd::Zero(height, width); + double tempvalue = 0; + for (int i = 0; i < height; i++) { + fin.read((char*)&tempvalue, sizeof(double)); this->AntPosition(i, 0) = tempvalue; //X + fin.read((char*)&tempvalue, sizeof(double)); this->AntPosition(i, 1) = tempvalue;//Y + fin.read((char*)&tempvalue, sizeof(double)); this->AntPosition(i, 2) = tempvalue;//Z + + // echo + for (int j = 0; j < width; j++) { + fin.read((char*)&tempvalue, sizeof(double)); // real + this->echo(i, j) = this->echo(i, j) + tempvalue; + } + //imag + for (int j = 0; j < width; j++) { + fin.read((char*)&tempvalue, sizeof(double)); // real + this->echo(i, j) = this->echo(i, j) + std::complex(0, tempvalue); + } + } + fin.close(); + return -1; +} + +int LAMPTOOLAPI BP2DProcessClass::saveTiFF(Eigen::MatrixXcd m) +{ + return saveMatrixXcd2TiFF(m, this->out_path); +} + +/// +/// 成像工作流，注意存在大量的内存浪费，后期可以根据情况进行优化 +/// +/// +int LAMPTOOLAPI BP2DProcessClass::start() +{ + // step 0 生成文件夹路径，为中间临时文件输出,构建临时环境，正式版需要注释相关代码 + QString parantPath = getParantFolderNameFromPath(this->out_path); + + + // step 1 插值频率域 -- FEKO 频率插值是均匀插值 + size_t Nf = this->freqnum; + size_t Nxa = this->AntPosition.rows(); + double df = this->Frequencylist(1) - this->Frequencylist(0); // FEKO 的频率插值比较规整 + double startfreq = this->Frequencylist(0); // 起始频率 + + // 输出频域 图像（未处理前） 并将复数数据转换为 振幅 与 相位 ，行：方位向，列：距离向------------------- + QString before_path = JoinPath(parantPath, "freq_ampDB_angle_echo.tiff"); + WriteComplexData2AmpdB_Arg(before_path, this->echo);// 保存为频域回波 + this->logFUN(100, "Read Echo\n"); + // ---------------------------------------------------------------------------------------------- + + // ---------------------------------------------------------------------------------------------- + + this->logFUN(100, "hamming windows over"); + // ----------------- step2 傅里叶变换 （频域-->时域） ----------------------------------------- + Eigen::MatrixXcd echoTime = IFFTW1D(this->echo); + size_t timefreqnum = echoTime.cols(); + double timeFreqBandWidth = df * (timefreqnum - 1); + this->logFUN(100, "echo IFFT over"); + + // 输出频域 图像（加窗） 并将复数数据转换为 振幅 与 相位 ，行：方位向，列：距离向------------------- + QString time_echo_path = JoinPath(parantPath, "freq_ampDB_angle_echo_hanning_ifft.tiff"); + WriteComplexData2AmpdB_Arg(time_echo_path, echoTime);// 保存为频域回波 + this->logFUN(100, "after ifft ,in Time\n"); + // ---------------------------------------------------------------------------------------------- + // step 3 时域域加窗 + Nf = echoTime.cols(); + Eigen::MatrixXd wfxa = Hanning(Nf, Nxa, 0.54); // wfxa = hanning(Nf,Nxa,alpha=0.54) + double normw = wfxa.array().sum(); // normw = TOTAL(wfxa) + + Eigen::MatrixXcd dfrex = (echoTime.array())* (wfxa.array()); // dfrex = dfrex * wfxa + // 输出频域 图像（加窗） 并将复数数据转换为 振幅 与 相位 ，行：方位向，列：距离向------------------- + QString hanning_path = JoinPath(parantPath, "freq_ampDB_angle_echo_hanning.tiff"); + WriteComplexData2AmpdB_Arg(hanning_path, dfrex);// 保存为频域回波 + this->logFUN(100, "after hanning\n"); + echoTime = dfrex; + // ----------------- step4 TBP成像 （时域） ----------------------------------------------------------- + Eigen::MatrixXcd img = BP2DImageByPixel(echoTime, this->AntPosition, this->minX, this->maxX, this->minY, this->maxY, this->Zplane, this->Rref, this->ImageWidth, this->ImageHeight, startfreq, timefreqnum, timeFreqBandWidth, this); + this->logFUN(100, "echo bp over"); + + + + // 输出 最终结果图像 并将复数数据转换为 振幅 与 相位 ，行：方位向，列：距离向------------------- + QString FBP_resultDB_path = JoinPath(parantPath, "freq_ampDB_angle_echo_hanning_BP.tiff"); + WriteComplexData2AmpdB_Arg(FBP_resultDB_path, img);// 保存为频域回波 + // ---------------------------------------------------------------------------------------- + // 输出 最终结果图像 并将复数数据转换为 振幅 与 相位 ，行：方位向，列：距离向------------------- + QString FBP_result_path = JoinPath(parantPath, "freq_amp_angle_echo_hanning_BP.tiff"); + WriteComplexData2Amp_Arg(FBP_result_path, img);// 保存为频域回波 + // ---------------------------------------------------------------------------------------- + + this->saveTiFF(img); + this->logFUN(100, "out bpimag over : " + this->out_path); + return -1; +} + +int LAMPTOOLAPI BP2DProcessClass::logFUN(int percent, QString logtext) { + qDebug() << "\rBPProcess [" << percent << "%]\t" << logtext; + if (percent < 100) { + qDebug()<<"\n"; + } + return 0; +}; + + + + + + +//////////////////////////////////////////////////////////////////////////////////////// +/// +/// FBP成像工作流（频域），参考 bp_linear_2d +/// Juan M Lopez-Sanchez, University of Alicante +/// +/// +//////////////////////////////////////////////////////////////////////////////////////// + + + + + + + + + + + + + + + + +int FBP2DProcessClass::initProcess(QString in_path, QString out_path, double Rref, double minX, double maxX, double minY, double maxY, double PlaneZ, int ImageHeight, int ImageWidth) +{ + // 初始化参数 + this->readEchoFile(in_path); + this->out_path = out_path; + this->Rref = Rref; // 成像中心的参考距离 + this->minX = minX; + this->maxX = maxX; + this->minY = minY; + this->maxY = maxY; + this->Zplane = PlaneZ; + this->ImageHeight = ImageHeight; + this->ImageWidth = ImageWidth; + // 计算坐标的到成像中心的坐标 + this->centerX = (this->maxX + this->minX) / 2; + this->centerY = (this->minY + this->maxY) / 2; + this->AntPosition.col(0) = this->AntPosition.col(0).array() - this->centerX; + this->AntPosition.col(1) = this->AntPosition.col(1).array() - this->centerY; + this->AntPosition.col(2) = this->AntPosition.col(2).array() - this->Zplane; + // 计算频率相关信息 + this->f1 = this->Frequencylist(0); + this->freqnum = this->Frequencylist.rows(); + this->fc = (this->f1 + this->Frequencylist(this->freqnum - 1)); + return -1; +} + +int FBP2DProcessClass::readEchoFile(QString in_path) +{ + std::ifstream fin(in_path.toUtf8().constData(), std::ios::in | std::ios::binary); + if (!fin.is_open()) exit(2); + // 读取参数 + int height = 0; + int width = 0; + fin.read((char*)&height, sizeof(int)); + fin.read((char*)&width, sizeof(int)); + this->height = height; + this->width = width; + //frequencylist[width*double] + double* freqs = new double[width]; + fin.read((char*)freqs, sizeof(double) * width); + this->Frequencylist = Eigen::VectorXd::Zero(width, 1); // 列向量 + for (int i = 0; i < width; i++) { + this->Frequencylist(i) = freqs[i]; + } + delete[] freqs; + freqs = nullptr; + + // 读取回波并解析脉冲天线位置 + this->AntPosition = Eigen::MatrixXd::Zero(height, 3); + this->echo = Eigen::MatrixXcd::Zero(height, width); + double tempvalue = 0; + for (int i = 0; i < height; i++) { + fin.read((char*)&tempvalue, sizeof(double)); this->AntPosition(i, 0) = tempvalue; //X + fin.read((char*)&tempvalue, sizeof(double)); this->AntPosition(i, 1) = tempvalue;//Y + fin.read((char*)&tempvalue, sizeof(double)); this->AntPosition(i, 2) = tempvalue;//Z + + // echo + for (int j = 0; j < width; j++) { + fin.read((char*)&tempvalue, sizeof(double)); // real + this->echo(i, j) = this->echo(i, j) + tempvalue; + } + //imag + for (int j = 0; j < width; j++) { + fin.read((char*)&tempvalue, sizeof(double)); // real + this->echo(i, j) = this->echo(i, j) + std::complex(0, tempvalue); + } + } + fin.close(); + return -1; +} + +int FBP2DProcessClass::saveTiFF(Eigen::MatrixXcd m) +{ + return saveMatrixXcd2TiFF(m, this->out_path); +} + + + + + + + + + + + + + + + + + +/// +/// FBP成像工作流（频域），参考 bp_linear_2d +/// Juan M Lopez-Sanchez, University of Alicante +/// +/// +int FBP2DProcessClass::start() +{ + // step 0 生成文件夹路径，为中间临时文件输出,构建临时环境，正式版需要注释相关代码 + QString parantPath = getParantFolderNameFromPath(this->out_path); + + + // step 1 插值频率域 -- FEKO 频率插值是均匀插值 + size_t Nf = this->freqnum; + size_t Nxa = this->AntPosition.rows(); + // 输出频域 图像（未处理前） 并将复数数据转换为 振幅 与 相位 ，行：方位向，列：距离向------------------- + QString before_path = JoinPath(parantPath, "freq_ampDB_angle_echo.tiff"); + WriteComplexData2AmpdB_Arg(before_path, this->echo);// 保存为频域回波 + this->logFUN(100, "Read Echo\n"); + // ---------------------------------------------------------------------------------------------- + // step 2 频率域加窗 + Eigen::MatrixXd wfxa = Hanning(Nf, Nxa, 0.54); // wfxa = hanning(Nf,Nxa,alpha=0.54) + double normw = wfxa.array().sum(); // normw = TOTAL(wfxa) + + Eigen::MatrixXcd dfrex = (this->echo.array()) * (wfxa.array()); // dfrex = dfrex * wfxa + // 输出频域 图像（加窗） 并将复数数据转换为 振幅 与 相位 ，行：方位向，列：距离向------------------- + QString hanning_path = JoinPath(parantPath, "freq_ampDB_angle_echo_hanning.tiff"); + WriteComplexData2AmpdB_Arg(hanning_path, dfrex);// 保存为频域回波 + this->logFUN(100, "after hanning\n"); + // ---------------------------------------------------------------------------------------------- + // step 3 Freq Backprojection algorithm + Eigen::MatrixXcd result_image = Eigen::MatrixXcd::Zero(this->ImageHeight, this->ImageWidth); // 构建网格，每个网格坐标根据中心点，平面点计算 + double dx = (this->maxX - this->minX) / (this->ImageWidth - 1); + double dy = (this->maxY - this->minY) / (this->ImageWidth - 1); + double px = 0, py = 0, pz = 0; // 像素坐标 + Eigen::MatrixXd Xdis = Eigen::MatrixXd::Zero(Nxa, Nf); + + // ; Nxa Nf + // mfre = fre # replicate(1, Nxa) ---- Nxa * Nf 复制每行 f1,f2,f3 ....... + // c = 0.2997925 ; Speed of light + // factorj = cj * (4.0*!PI/c) * mfre --- Nxa * Nf 复制每行 f1,f2,f3 ....... + // + + Eigen::MatrixXcd factorj = Eigen::MatrixXcd::Zero(Nxa, Nf); + for (size_t i = 0; i < Nxa; i++) { + for (size_t j = 0; j < Nf; j++) { + factorj(i, j) = std::complex(0, 4.0 * PI / (LIGHTSPEED * 1e-9)) * (this->Frequencylist(j) * 1e-9); // cj * (4.0*!PI/c) * mfre + } + } + + double R0 = this->Rref;// -- 参考斜距 + Eigen::MatrixXd R = Eigen::MatrixXd::Zero(Nxa, Nf); + Eigen::MatrixXd distR = Eigen::MatrixXd::Zero(Nxa, 1); + Eigen::MatrixXcd term = Eigen::MatrixXcd::Zero(Nxa, Nf); + Eigen::MatrixXcd data = Eigen::MatrixXcd::Zero(Nxa, Nf); + + for (size_t ix = 0; ix < this->ImageWidth; ix++) { + if (ix % 100 == 0) { + this->logFUN(size_t(ix * 1.0 / this->ImageWidth * 100), "FBP ....\n"); + } + // theta = theta / !radeg -- 根据入射角 计算 参考坐标 + // ya = Ro * sin(theta) -- 雷达 y 坐标 + // za = Ro * cos(theta) -- 雷达 z 坐标 + // + // xa_min = pos_start + // xa_max = pos_stop + // Nxa = Npos + // + // dxa = ( xa_max - xa_min ) / float(Nxa-1) + // xa = xa_min + findgen(Nxa) * dxa -- 雷达 x 坐标 + // + // dx = (x_max-x_min) / float(Nx-1) + // x = x_min + findgen(Nx) * dx + // + // dy = (y_max-y_min) / float(Ny-1) + // y = y_min + findgen(Ny) * dy + // mxa = replicate(1, Nf) # xa ;--- Nxa * Nf 每个列 ant_x1， ant_x2, ant_x3,...... + // xdis = (x[ix]-mxa)^2 + za^2 ;--- Nxa * Nf 每个列 ant_x1_z1，ant_x2_z2,ant_x3_z3,...... + // R = sqrt( xdis + (ya - y[iy])^2 ) ;--- Nxa * Nf 每个列 ant_x1， ant_x2, ant_x3,...... + // + px = ix * dx + this->minX; + pz = this->Zplane; + for (size_t iy = 0; iy < this->ImageHeight; iy++) { + py = iy * dy + this->minY; + + // mxa = replicate(1, Nf) # xa ;--- Nxa * Nf 每个列 ant_x1， ant_x2, ant_x3,...... + // xdis = (x[ix]-mxa)^2 + za^2 ;--- Nxa * Nf 每个列 ant_x1_z1，ant_x2_z2,ant_x3_z3,...... + // R = sqrt( xdis + (ya - y[iy])^2 ) ;--- Nxa * Nf 每个列 ant_x1， ant_x2, ant_x3,...... + // + + distR = ((this->AntPosition.col(0).array() - px).pow(2) + (this->AntPosition.col(1).array() - py).pow(2) + (this->AntPosition.col(2).array() - pz).pow(2)).array().sqrt().array(); // 复制 + + for (size_t jj = 0; jj < Nf; jj++) { // 逐列复制 + R.col(jj) = distR.array(); + } + // term = (R/Ro)^2 * exp( factorj * ( R - Ro ) ) ; Nxa x Nf + + term = (R.array() / R0).array().pow(2) * ((factorj.array() * (R.array() - R0)).array().exp()); + data = dfrex.array() * term.array(); + result_image(iy, ix) = data.array().sum(); + // 临时文件输出----------------------- + if (iy == this->ImageHeight / 2 && ix == this->ImageWidth / 2) { + // 输出 校正值 并将复数数据转换为 振幅 与 相位 ，行：方位向，列：距离向------------------- + QString FBP_termDB_path = JoinPath(parantPath, "freq_ampDB_angle_echo_hanning_term_FBP.tiff"); + WriteComplexData2AmpdB_Arg(FBP_termDB_path, term);// 保存为频域回波 + + // 输出 校正后图像 并将复数数据转换为 振幅 与 相位 ，行：方位向，列：距离向------------------- + QString FBP_dataDB_path = JoinPath(parantPath, "freq_ampDB_angle_echo_hanning_data_FBP.tiff"); + WriteComplexData2AmpdB_Arg(FBP_dataDB_path, data);// 保存为频域回波 + + // ---------------------------------------------------------------------------------------------- + } + + /***/ + } + } + result_image = result_image.array() / normw; + // 输出 最终结果图像 并将复数数据转换为 振幅 与 相位 ，行：方位向，列：距离向------------------- + QString FBP_resultDB_path = JoinPath(parantPath, "freq_ampDB_angle_echo_hanning_FBP.tiff"); + WriteComplexData2AmpdB_Arg(FBP_resultDB_path, result_image);// 保存为频域回波 + // ---------------------------------------------------------------------------------------------- + // 输出 最终结果图像 并将复数数据转换为 振幅 与 相位 ，行：方位向，列：距离向------------------- + QString FBP_result_path = JoinPath(parantPath, "freq_amp_angle_echo_hanning_FBP.tiff"); + WriteComplexData2Amp_Arg(FBP_result_path, result_image);// 保存为频域回波 + // ---------------------------------------------------------------------------------------------- + this->saveTiFF(result_image); + + return 0; +} + +int FBP2DProcessClass::logFUN(int percent, QString logtext) +{ + qDebug() << "\rFBPProcess [" << percent << "%]\t" << logtext; + if (percent < 100) { + qDebug() << "\n"; + } + return 0; +} + +int FEKOFarFieldProcessClass::initProcess(QString in_path, QString out_path, double Rref, double minX, double maxX, double minY, double maxY, double PlaneZ, int ImageHeight, int ImageWidth) +{ + // 初始化参数 + this->readEchoFile(in_path); + this->out_path = out_path; + this->Rref = Rref; // 成像中心的参考距离 + this->minX = minX; + this->maxX = maxX; + this->minY = minY; + this->maxY = maxY; + this->Zplane = PlaneZ; + this->ImageHeight = ImageHeight; + this->ImageWidth = ImageWidth; + // 计算坐标的到成像中心的坐标 + this->centerX = (this->maxX + this->minX) / 2; + this->centerY = (this->minY + this->maxY) / 2; + this->AntPosition.col(0) = this->AntPosition.col(0).array() - this->centerX; + this->AntPosition.col(1) = this->AntPosition.col(1).array() - this->centerY; + this->AntPosition.col(2) = this->AntPosition.col(2).array() - this->Zplane; + // 计算频率相关信息 + this->f1 = this->Frequencylist(0); + this->freqnum = this->Frequencylist.rows(); + this->fc = (this->f1 + this->Frequencylist(this->freqnum - 1)); + return -1; +} + +int FEKOFarFieldProcessClass::readEchoFile(QString in_path) +{ + std::ifstream fin(in_path.toUtf8().constData(), std::ios::in | std::ios::binary); + if (!fin.is_open()) exit(2); + // 读取参数 + int height = 0; + int width = 0; + fin.read((char*)&height, sizeof(int)); + fin.read((char*)&width, sizeof(int)); + this->height = height; + this->width = width; + //frequencylist[width*double] + double* freqs = new double[width]; + fin.read((char*)freqs, sizeof(double) * width); + this->Frequencylist = Eigen::VectorXd::Zero(width, 1); // 列向量 + for (int i = 0; i < width; i++) { + this->Frequencylist(i) = freqs[i]; + } + delete[] freqs; + freqs = nullptr; + + // 读取回波并解析脉冲天线位置 + this->AntPosition = Eigen::MatrixXd::Zero(height, 3); + this->echo = Eigen::MatrixXcd::Zero(height, width); + double tempvalue = 0; + for (int i = 0; i < height; i++) { + fin.read((char*)&tempvalue, sizeof(double)); this->AntPosition(i, 0) = tempvalue; //X + fin.read((char*)&tempvalue, sizeof(double)); this->AntPosition(i, 1) = tempvalue;//Y + fin.read((char*)&tempvalue, sizeof(double)); this->AntPosition(i, 2) = tempvalue;//Z + + // echo + for (int j = 0; j < width; j++) { + fin.read((char*)&tempvalue, sizeof(double)); // real + this->echo(i, j) = this->echo(i, j) + tempvalue; + } + //imag + for (int j = 0; j < width; j++) { + fin.read((char*)&tempvalue, sizeof(double)); // real + this->echo(i, j) = this->echo(i, j) + std::complex(0, tempvalue); + } + } + fin.close(); + return -1; +} + +int FEKOFarFieldProcessClass::saveTiFF(Eigen::MatrixXcd m) +{ + return saveMatrixXcd2TiFF(m, this->out_path); +} + +int FEKOFarFieldProcessClass::start() +{ + + //// step 0 生成文件夹路径，为中间临时文件输出,构建临时环境，正式版需要注释相关代码 + //QString parantPath = getParantFolderNameFromPath(this->out_path); + + + //double freqbandWidth = this->Frequencylist(this->get_Nf() - 1) - this->Frequencylist(0); + //double dx = LIGHTSPEED / 2.0 / std::abs(freqbandWidth); // 距离向分辨率 + //double startPhiAngle = this->AntPosition(0, 1); + //double endPhiAngle = this->AntPosition(this->get_Nxa(), 1); + //double AzAngleWidth = std::abs(endPhiAngle - startPhiAngle); + //double halfAzAngleWidth = AzAngleWidth / 2.0; + //double AzBandWidth = 2 * std::tan(Degrees2Radians(halfAzAngleWidth)) * this->get_minFreq(); + //double dy = LIGHTSPEED / 2.0 / std::abs(AzBandWidth); // 方位向分辨率 + + + + //// step 1 插值频率域 -- FEKO 频率插值是均匀插值 + //size_t Nf = this->freqnum; + //size_t Nxa = this->AntPosition.rows(); + + + + //double angStart = startPhiAngle; + //double angEnd = endPhiAngle; + //double freqEnd = this->Frequencylist(this->get_Nf() - 1); + //double fxStart = this->Frequencylist(0); + //double fyStart = std::tan(-1 * std::abs(halfAzAngleWidth)) * fxStart; + //double fyEnd = std::tan(-1 * std::abs(halfAzAngleWidth)) * fxStart; + //double fxEnd = std::sqrt(std::pow(freqEnd, 2) - std::pow(fyStart, 2)); + + //double freqSamples = this->Frequencylist.rows(); + //double angSamples = this->AntPosition.rows(); + + //double ysd = (fyEnd - fyStart) / (angSamples - 1); // 转换为 网格节点 + //double xsd = (fxEnd - fxStart) / (freqSamples - 1); // 转换为 网格节点 + + //Eigen::MatrixXd image_points_freq(angSamples, freqSamples); + //Eigen::MatrixXd image_points_ang(angSamples, freqSamples); + //Eigen::MatrixXcd Echo_data(angSamples, freqSamples); + + //Eigen::MatrixXd xi(Nxa,Nf); + //Eigen::MatrixXd yi(Nxa, Nf); + + + //// 创建网格 + //for (int f = 0; f < freqSamples; f++) { + // for (int a = 0; a < angSamples; a++) { + // double x = fxStart + f * xsd; + // double y = fyStart + a * ysd; + // image_points_freq(a, f) = std::sqrt(x * x + y * y); // 获取每个点的 频率 + // image_points_ang(a, f) = std::atan2(y, x); // 获取每个点的 角度坐标 + // Echo_data(a, f) = InterpolateComplex(xi, yi, this->echo, image_points_freq(a, f), image_points_ang(a, f)); // 插值得到区域数据 + + // } + //} + // + //// step 2 频率域加窗 + //Eigen::MatrixXd wfxa = Hanning(Nf, Nxa, 0.54); // wfxa = hanning(Nf,Nxa,alpha=0.54) + //double normw = wfxa.array().sum(); // normw = TOTAL(wfxa) + + //Eigen::MatrixXcd dfrex = (this->echo.array()) * (wfxa.array()); // dfrex = dfrex * wfxa + //// 输出频域 图像（加窗） 并将复数数据转换为 振幅 与 相位 ，行：方位向，列：距离向------------------- + //QString hanning_path = JoinPath(parantPath, "freq_ampDB_angle_echo_hanning.tiff"); + //WriteComplexData2AmpdB_Arg(hanning_path, dfrex);// 保存为频域回波 + //this->logFUN(100, "after hanning\n"); + + //Echo_data = dfrex; + + // + //// 执行二维傅立叶变换 + //Eigen::MatrixXcd H_echo_win_fft = FFTW2D(Echo_data); + + //// 获取矩阵的行数和列数 + //int n_freq = Echo_data.rows(); + //int n_angle = Echo_data.cols(); + + //// 将傅立叶变换结果归一化 + //H_echo_win_fft = H_echo_win_fft.array()/(n_freq * n_angle); + + //// 执行逆傅立叶变换（如果需要） + //// MatrixXcd H_echo_win_ifft = ifft2(H_echo_win_fft); + //// H_echo_win_ifft /= (n_freq * n_angle); + + //// 执行平移 + //Eigen::MatrixXcd result_image = fftshift(H_echo_win_fft); + + + //// 输出 最终结果图像 并将复数数据转换为 振幅 与 相位 ，行：方位向，列：距离向------------------- + //QString FBP_resultDB_path = JoinPath(parantPath, "freq_ampDB_angle_echo_hanning_FBP.tiff"); + //WriteComplexData2AmpdB_Arg(FBP_resultDB_path, result_image);// 保存为频域回波 + //// ---------------------------------------------------------------------------------------------- + //// 输出 最终结果图像 并将复数数据转换为 振幅 与 相位 ，行：方位向，列：距离向------------------- + //QString FBP_result_path = JoinPath(parantPath, "freq_amp_angle_echo_hanning_FBP.tiff"); + //WriteComplexData2Amp_Arg(FBP_result_path, result_image);// 保存为频域回波 + //// ---------------------------------------------------------------------------------------------- + //this->saveTiFF(result_image); + + + return 0; +} + +int FEKOFarFieldProcessClass::logFUN(int percent, QString logtext) +{ + qDebug() << "\rFBPProcess [" << percent << "%]\t" << logtext; + if (percent < 100) { + qDebug() << "\n"; + } + return 0; +} diff --git a/src/LAMPTool/FEKONearBPBasic.h b/src/LAMPTool/FEKONearBPBasic.h new file mode 100644 index 0000000..8a634fd --- /dev/null +++ b/src/LAMPTool/FEKONearBPBasic.h @@ -0,0 +1,161 @@ +#pragma once +#ifndef _FEKONEARBPBASIC_H_ +#define _FEKONEARBPBASIC_H_ + +/** +* 适用于FEKO的近场结果的 BP 成像算法 +* +**/ +#include "LAMPToolAPI.h" +#include +#include +#include +#include +#include +#include +#include +#include "referenceHeader.h" + + + +// FEKO 几何关系处理 + +/// +/// 笛卡尔坐标系，转换为 极坐标系 +/// +/// [X,Y,Z;X1,Y1,Z1] +/// +Eigen::MatrixXd LAMPTOOLAPI Cartesian2Spherical(Eigen::MatrixXd CartesianPoint); + +/// +/// 时域BP +/// +class LAMPTOOLAPI BP2DProcessClass { +public: + size_t height; + size_t width; + Eigen::MatrixXcd echo; + Eigen::MatrixXd AntPosition; + Eigen::VectorXd Frequencylist; + size_t freqnum; + double f1; + double fc; + double Rref; // 成像中心的参考距离 + double minX; + double maxX; + double minY; + double maxY; + double centerX; + double centerY; + double Zplane; + size_t ImageHeight; + size_t ImageWidth; + QString out_path; +public: + virtual int initProcess(QString in_path, QString out_path, double Rref, double minX, double maxX, double minY, double maxY, double PlaneZ, int ImageHeight, int ImageWidth); + virtual int readEchoFile(QString in_path); + virtual int saveTiFF(Eigen::MatrixXcd m); + virtual int start(); + virtual int logFUN(int percent, QString logtext); +}; + + +int LAMPTOOLAPI BP2DProcess(QString in_path, QString out_path, double Rref, double minX, double maxX, double minY, double maxY, double PlaneZ,int ImageHeight,int ImageWidth); +// BP 成像时，逐像素点计算，计算速度慢，回波（PRFNUM,freqNUM) +template +Eigen::MatrixXcd LAMPTOOLAPI BP2DImageByPixel(Eigen::MatrixXcd timeEcho,Eigen::MatrixXd AntPosition,double minX,double maxX,double minY,double maxY,double PlaneZ ,double Rref,size_t ImageWidth,size_t ImageHeight,double startfreq,size_t timefreqnum,double timeFreqBandWidth, T* logclss=nullptr); + +/// +/// 时域BP --- FBP +/// +class LAMPTOOLAPI FBP2DProcessClass:public BP2DProcessClass { +public: + size_t height; + size_t width; + Eigen::MatrixXcd echo; + Eigen::MatrixXd AntPosition; + Eigen::VectorXd Frequencylist; + size_t freqnum; + double f1; + double fc; + double Rref; // 成像中心的参考距离 + double minX; + double maxX; + double minY; + double maxY; + double centerX; + double centerY; + double Zplane; + size_t ImageHeight; + size_t ImageWidth; + QString out_path; +public: + virtual int initProcess(QString in_path, QString out_path, double Rref, double minX, double maxX, double minY, double maxY, double PlaneZ, int ImageHeight, int ImageWidth); + virtual int readEchoFile(QString in_path); + virtual int saveTiFF(Eigen::MatrixXcd m); + virtual int start(); + virtual int logFUN(int percent, QString logtext); +}; + +// BP 成像时，逐脉冲计算，回波（PRFNUM,freqNUM) +template +Eigen::MatrixXcd LAMPTOOLAPI BP2DImageByPluse(Eigen::MatrixXcd timeEcho, Eigen::MatrixXd AntPosition, double minX, double maxX, double minY, double maxY, double PlaneZ, double Rref, size_t ImageWidth, size_t ImageHeight, double startfreq, size_t timefreqnum, double timeFreqBandWidth, T* logclss = nullptr); + +int LAMPTOOLAPI FBP2DProcess(QString in_path, QString out_path, double Rref, double minX, double maxX, double minY, double maxY, double PlaneZ, int ImageHeight, int ImageWidth); + +// 生成成像的数据文件 +int LAMPTOOLAPI build2Bin(QString path,int width,int height,std::vector &freqs,Eigen::MatrixXd AntPositions,Eigen::MatrixXcd echo); + + +// 定义插值函数，以处理复数值 +std::complex LAMPTOOLAPI InterpolateComplex(Eigen::MatrixXd& xi, Eigen::MatrixXd& yi, Eigen::MatrixXcd& data, double x, double y); + + + +/// +/// 远场成像 +/// +class LAMPTOOLAPI FEKOFarFieldProcessClass :public BP2DProcessClass { +public: + size_t height; + size_t width; + Eigen::MatrixXcd echo; + Eigen::MatrixXd AntPosition; // theta phi 0 + Eigen::VectorXd Frequencylist; + size_t freqnum; + double f1; + double fc; + double Rref; // 成像中心的参考距离 + double minX; + double maxX; + double minY; + double maxY; + double centerX; + double centerY; + double Zplane; + size_t ImageHeight; + size_t ImageWidth; + QString out_path; +public: + virtual int initProcess(QString in_path, QString out_path, double Rref, double minX, double maxX, double minY, double maxY, double PlaneZ, int ImageHeight, int ImageWidth); + virtual int readEchoFile(QString in_path); + virtual int saveTiFF(Eigen::MatrixXcd m); + virtual int start(); + virtual int logFUN(int percent, QString logtext); + + // 平均频率间隔 + double get_df() { return (this->Frequencylist(this->get_Nf() - 1) - this->Frequencylist(0)) / (this->get_Nf() - 1); }; + // 脉冲数 + size_t get_Nxa() { return this->AntPosition.rows(); }; + // 频率点数 + size_t get_Nf() { return this->Frequencylist.rows(); }; + double get_minFreq() { return this->Frequencylist(this->get_Nf() - 1) > this->Frequencylist(0) ? this->Frequencylist(0) : this->Frequencylist(this->get_Nf() - 1); }; + +}; + + + +#endif + + + diff --git a/src/LAMPTool/FEKONearBpBaseImage.cpp b/src/LAMPTool/FEKONearBpBaseImage.cpp new file mode 100644 index 0000000..ff3e27f --- /dev/null +++ b/src/LAMPTool/FEKONearBpBaseImage.cpp @@ -0,0 +1,100 @@ +// FEKONearBpBaseImage.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。 +// +#include "FEKONearBPBasic.h" +#include +#include +#include "LampToolTest.h" + +int TestImageBP_main(int argc, char* argv[]) +{ + qDebug() << "process starting .......\n"; + qDebug() << "-i in_echo.bin -o out_resule.bin -Rref 10 -minX -1 -maxX 1 -minY -1 -maxY 1 -PlaneZ 0 -ImageHeight 201 -ImageWidth 201" << "\n"; + qDebug() << "BP mehtod for FEKO Near field (Cartesian) " << "\n"; + qDebug() << "format: " << "\n"; + qDebug() << "height[int32] width[int32] frequencylist[width*double] echoData " << "\n"; + qDebug() << "echoData format: ant[x,y,z 3xdouble] real[width*double] imag[width*double]" << "\n"; + qDebug() << "echodata [freqNUM,PRF_num]" << "\n"; + qDebug() << "warnning:Eigen default matrix is row-major storage" << "\n"; + qDebug() << "===============params list==============================" << "\n"; + QString in_path = ""; + QString out_path = ""; + int mode = 1; + double Rref=0, minX=0, maxX = 0, minY = 0, maxY = 0, PlaneZ = 0, ImageHeight = 0, ImageWidth = 0; + for (int i = 0; i < argc; i++) { + qDebug() << argv[i] << "\n"; + } + for (int i = 0; i < argc; i++) { + if (strcmp(argv[i], "-i")==0) { + if (i == argc - 1) { return 1; } + in_path = argv[i + 1]; + qDebug() << argv[i] << " = " << argv[i + 1] << "\n"; + } + else if (strcmp(argv[i], "-o") == 0) { + if (i == argc - 1) { return 1; } + out_path = argv[i + 1]; + qDebug() << argv[i] << " = " << argv[i + 1] << "\n"; + } + else if (strcmp(argv[i], "-Rref") == 0) { + if (i == argc - 1) { return 1; } + Rref =std::stod(argv[i+1]); + qDebug() << argv[i] << " = " << argv[i + 1] << "\n"; + } + else if (strcmp(argv[i], "-minX") == 0) { + if (i == argc - 1) { return 1; } + minX = std::stod(argv[i + 1]); + qDebug() << argv[i] << " = " << argv[i + 1] << "\n"; + } + else if (strcmp(argv[i], "-maxX") == 0) { + if (i == argc - 1) { return 1; } + maxX = std::stod(argv[i + 1]); + qDebug() << argv[i] << " = " << argv[i + 1] << "\n"; + } + else if (strcmp(argv[i], "-minY") == 0) { + if (i == argc - 1) { return 1; } + minY = std::stod(argv[i + 1]); + qDebug() << argv[i] << " = " << argv[i + 1] << "\n"; + } + else if (strcmp(argv[i], "-maxY") == 0) { + if (i == argc - 1) { return 1; } + maxY = std::stod(argv[i + 1]); + qDebug() << argv[i] << " = " << argv[i + 1] << "\n"; + } + else if (strcmp(argv[i], "-PlaneZ") == 0) { + if (i == argc - 1) { return 1; } + PlaneZ = std::stod(argv[i + 1]); + qDebug() << argv[i] << " = " << argv[i + 1] << "\n"; + } + else if (strcmp(argv[i], "-ImageHeight") == 0) { + if (i == argc - 1) { return 1; } + ImageHeight = std::stod(argv[i + 1]); + qDebug() << argv[i] << " = " << argv[i + 1] << "\n"; + } + else if (strcmp(argv[i], "-ImageWidth")==0) { + if (i == argc - 1) { return 1; } + ImageWidth = std::stod(argv[i + 1]); + qDebug() << argv[i] << " = " << argv[i + 1] << "\n"; + } + else if (strcmp(argv[i], "-mode") == 0) { + if (i == argc - 1) { return 1; } + mode = std::stoi(argv[i + 1]); + qDebug() << argv[i] << " = " << argv[i + 1] << "\n"; + } + } + qDebug() << "======================================================" << "\n"; + //qDebug() << Eigen::MatrixXd::LinSpaced(1, 10).array() << "\n"; + if (mode == 0) { + //BP2DProcess(in_path, out_path, Rref, minX, maxX, minY, maxY, PlaneZ, ImageHeight, ImageWidth); + } + else if (mode == 1) { + //FBP2DProcess(in_path, out_path, Rref, minX, maxX, minY, maxY, PlaneZ, ImageHeight, ImageWidth); + } + else if (mode == 2) { + + } + + return 0; +} + + + + diff --git a/src/LAMPTool/FileOperator.cpp b/src/LAMPTool/FileOperator.cpp new file mode 100644 index 0000000..8aac6f7 --- /dev/null +++ b/src/LAMPTool/FileOperator.cpp @@ -0,0 +1,188 @@ +#include "FileOperator.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +std::vector LAMPTOOLAPI getFilelist(const QString& folderpath, const QString& filenameExtension, int (*logfun)(QString logtext, int value)) +{ + QString filenameExtensionStr = filenameExtension; + filenameExtensionStr = filenameExtensionStr.remove(0, 1); + std::vector filenames(0); + if (isExists(folderpath) && isDirectory(folderpath)) { + QDir directory(folderpath); + if (directory.exists() && directory.isReadable()) { + QFileInfoList fileList = directory.entryInfoList(QDir::Files | QDir::NoDotAndDotDot); + for (const QFileInfo& fileInfo : fileList) { + qDebug() << fileInfo.filePath() + "\nExtension: (" + filenameExtensionStr + ", " + fileInfo.suffix() + ")"; + if (filenameExtensionStr == u8"*" || filenameExtensionStr == fileInfo.suffix()) { + filenames.push_back(fileInfo.filePath()); + } + if (logfun) { + logfun(fileInfo.filePath() + "\nExtension: (" + filenameExtensionStr + ", " + fileInfo.suffix() + ")", 4); + } + } + } + else { + qWarning() << "Folder does not exist or is not readable: " << folderpath; + } + return filenames; + } + else { + return std::vector(0); + } + +} + +QString LAMPTOOLAPI getParantFolderNameFromPath(const QString& path) +{ + QDir directory(path); + directory.cdUp(); + QString parentPath = directory.absolutePath(); + return directory.dirName(); +} + +QString LAMPTOOLAPI getParantFromPath(const QString& path) +{ + //qDebug() << path; + QDir directory(path); + directory.cdUp(); + QString parentPath = directory.absolutePath(); + //qDebug() << parentPath; + return parentPath; +} + +QString LAMPTOOLAPI getFileNameFromPath(const QString& path) +{ + QFileInfo fileInfo(path); + return fileInfo.fileName(); +} + +bool LAMPTOOLAPI isDirectory(const QString& path) +{ + QFileInfo fileinfo(path); + return fileinfo.isDir(); +} + +bool LAMPTOOLAPI isExists(const QString& path) +{ + QFileInfo fileinfo(path); + return fileinfo.exists(); +} + +bool LAMPTOOLAPI isFile(const QString& path) +{ + QFileInfo fileinfo(path); + return fileinfo.isFile(); +} + +int LAMPTOOLAPI write_binfile(char* filepath, char* data, size_t data_len) +{ + FILE* pd = fopen(filepath, "w"); + if (NULL == pd) { + return 2; + } + //数据块首地址: "&a"，元素大小: "sizeof(unsigned __int8)"， 元素个数: "10"， 文件指针："pd" + fwrite(data, sizeof(char), data_len, pd); + fclose(pd); + return -1; +} + +LAMPTOOLAPI char* read_textfile(char* text_path, int* length) +{ + char* data = NULL; + FILE* fp1 = fopen(text_path, "r"); + if (fp1 == NULL) { + return NULL; + } + else {} + // 读取文件 + fseek(fp1, 0, SEEK_END); + int data_length = ftell(fp1); + data = (char*)malloc((data_length + 1) * sizeof(char)); + rewind(fp1); + if (data_length == fread(data, sizeof(char), data_length, fp1)) { + data[data_length] = '\0'; // 文件尾 + } + else { + free(data); + fclose(fp1); + return NULL; + } + fclose(fp1); + *length = data_length + 1; + return data; +} + +bool LAMPTOOLAPI exists_test(const QString& name) +{ + return isExists(name); +} + +size_t LAMPTOOLAPI fsize(FILE* fp) +{ + size_t n; + fpos_t fpos; // 当前位置 + fgetpos(fp, &fpos); // 获取当前位置 + fseek(fp, 0, SEEK_END); + n = ftell(fp); + fsetpos(fp, &fpos); // 恢复之前的位置 + return n; +} + + +void LAMPTOOLAPI removeFile(const QString& filePath) +{ + QFile file(filePath); + + if (file.exists()) { + if (file.remove()) { + qDebug() << "File removed successfully: " << filePath; + } + else { + qWarning() << "Failed to remove file: " << filePath; + } + } + else { + qDebug() << "File does not exist: " << filePath; + } +} + +unsigned long LAMPTOOLAPI convertToULong(const QString& input) { + bool ok; // Used to check if the conversion was successful + unsigned long result = input.toULong(&ok); + + if (!ok) { + qWarning() << "Conversion to unsigned long failed for input: " << input; + } + + return result; +} + + + +void LAMPTOOLAPI copyFile(const QString& sourcePath, const QString& destinationPath) { + QFile sourceFile(sourcePath); + QFile destinationFile(destinationPath); + + if (sourceFile.exists()) { + if (sourceFile.copy(destinationPath)) { + // 复制成功 + //QMessageBox::information(nullptr, u8"成功", u8"文件复制成功"); + } + else { + // 复制失败 + QMessageBox::critical(nullptr, QObject::tr("error"), QObject::tr("file copy error")); + } + } + else { + // 源文件不存在 + QMessageBox::warning(nullptr, QObject::tr("warning"), QObject::tr("Source file not found")); + } +} \ No newline at end of file diff --git a/src/LAMPTool/FileOperator.h b/src/LAMPTool/FileOperator.h new file mode 100644 index 0000000..ed76fd1 --- /dev/null +++ b/src/LAMPTool/FileOperator.h @@ -0,0 +1,50 @@ +#pragma once + +#ifndef FILEOPERATOR_H +#define FILEOPERATOR_H + +#include "../LAMPTool/referenceHeader.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +bool LAMPTOOLAPI isDirectory(const QString& path); +bool LAMPTOOLAPI isExists(const QString& path); +bool LAMPTOOLAPI isFile(const QString& path); +void LAMPTOOLAPI removeFile(const QString& filePath); +unsigned long LAMPTOOLAPI convertToULong(const QString& input); +/// +/// 获取文件(绝对路径） +/// +/// +/// +/// +std::vector LAMPTOOLAPI getFilelist(const QString& folderpath, const QString& FilenameExtension = ".*",int (*logfun)(QString logtext,int value)=nullptr); + +QString LAMPTOOLAPI getParantFolderNameFromPath(const QString& path); + +QString LAMPTOOLAPI getFileNameFromPath(const QString& path); + +int LAMPTOOLAPI write_binfile(char* filepath, char* data, size_t data_len); + +LAMPTOOLAPI char* read_textfile(char* text_path, int* length); + +bool LAMPTOOLAPI exists_test(const QString& name); + +size_t LAMPTOOLAPI fsize(FILE* fp); + +QString LAMPTOOLAPI getParantFromPath(const QString& path); +void LAMPTOOLAPI copyFile(const QString& sourcePath, const QString& destinationPath); +// QT FileOperator +#endif \ No newline at end of file diff --git a/src/LAMPTool/GeoOperator.cpp b/src/LAMPTool/GeoOperator.cpp new file mode 100644 index 0000000..fdc77ce --- /dev/null +++ b/src/LAMPTool/GeoOperator.cpp @@ -0,0 +1,272 @@ +#include "GeoOperator.h" +#include +#include +#include +#include +//#include +#include +#include +#include +#include +#include +#include +#include +#include +//#include +#include //#include "ogrsf_frmts.h" +#include +#include +#include "GeoOperator.h" + +using namespace std; +using namespace Eigen; + + +Landpoint LAMPTOOLAPI operator +(const Landpoint& p1, const Landpoint& p2) +{ + return Landpoint{ p1.lon + p2.lon,p1.lat + p2.lat,p1.ati + p2.ati }; +} + +Landpoint LAMPTOOLAPI operator -(const Landpoint& p1, const Landpoint& p2) +{ + return Landpoint{ p1.lon - p2.lon,p1.lat - p2.lat,p1.ati - p2.ati }; +} + +bool LAMPTOOLAPI operator ==(const Landpoint& p1, const Landpoint& p2) +{ + return p1.lat == p2.lat && p1.lon == p2.lon && p1.ati == p2.ati; +} + + + +Landpoint LAMPTOOLAPI operator *(const Landpoint& p, double scale) +{ + return Landpoint{ + p.lon * scale, + p.lat * scale, + p.ati * scale + }; +} + + +Landpoint LAMPTOOLAPI LLA2XYZ(const Landpoint& LLA) { + double L = LLA.lon * d2r; + double B = LLA.lat * d2r; + double H = LLA.ati; + + double sinB = sin(B); + double cosB = cos(B); + + //double N = a / sqrt(1 - e * e * sin(B) * sin(B)); + double N = a / sqrt(1 - eSquare * sinB * sinB); + Landpoint result = { 0,0,0 }; + result.lon = (N + H) * cosB * cos(L); + result.lat = (N + H) * cosB * sin(L); + //result.z = (N * (1 - e * e) + H) * sin(B); + result.ati = (N * (1 - 1 / f_inverse) * (1 - 1 / f_inverse) + H) * sinB; + return result; +} + + +Eigen::MatrixXd LAMPTOOLAPI LLA2XYZ(Eigen::MatrixXd landpoint) +{ + landpoint.col(0) = landpoint.col(0).array() * d2r; // lon + landpoint.col(1) = landpoint.col(1).array() * d2r; // lat + + Eigen::MatrixXd sinB = (landpoint.col(1).array().sin());//lat + Eigen::MatrixXd cosB = (landpoint.col(1).array().cos());//lat + + Eigen::MatrixXd N = a / ((1 - sinB.array().pow(2) * eSquare).array().sqrt()); + Eigen::MatrixXd result(landpoint.rows(), 3); + + result.col(0) = (N.array() + landpoint.col(2).array()) * cosB.array() * Eigen::cos(landpoint.col(0).array()).array(); //x + result.col(1) = (N.array() + landpoint.col(2).array()) * cosB.array() * Eigen::sin(landpoint.col(0).array()).array(); //y + + result.col(2) = (N.array() * (1 - 1 / f_inverse) * (1 - 1 / f_inverse) + landpoint.col(2).array()) * sinB.array(); //z + + return result; +} + + +Landpoint LAMPTOOLAPI XYZ2LLA(const Landpoint& XYZ) { + double tmpX = XYZ.lon;// + double temY = XYZ.lat;// + double temZ = XYZ.ati; + + double curB = 0; + double N = 0; + double sqrtTempXY = sqrt(tmpX * tmpX + temY * temY); + double calB = atan2(temZ, sqrtTempXY); + int counter = 0; + double sinCurB = 0; + while (abs(curB - calB) * r2d > epsilon && counter < 25) + { + curB = calB; + sinCurB = sin(curB); + N = a / sqrt(1 - eSquare * sinCurB * sinCurB); + calB = atan2(temZ + N * eSquare * sinCurB, sqrtTempXY); + counter++; + } + + Landpoint result = { 0,0,0 }; + result.lon = atan2(temY, tmpX) * r2d; + result.lat = curB * r2d; + result.ati = temZ / sinCurB - N * (1 - eSquare); + return result; +} + + + + +double LAMPTOOLAPI getAngle(const Landpoint& a, const Landpoint& b) +{ + double c = dot(a, b) / (getlength(a) * getlength(b)); + if (a.lon * b.lat - a.lat * b.lon >= 0) { + return acos(c > 1 ? 1 : c < -1 ? -1 : c) * r2d; + } + else { + return 360 - acos(c > 1 ? 1 : c < -1 ? -1 : c) * r2d; + } +} + +double LAMPTOOLAPI dot(const Landpoint& p1, const Landpoint& p2) +{ + return p1.lat * p2.lat + p1.lon * p2.lon + p1.ati * p2.ati; +} + +double LAMPTOOLAPI getlength(const Landpoint& p1) { + + return sqrt(dot(p1, p1)); +} + +Landpoint LAMPTOOLAPI crossProduct(const Landpoint& a, const Landpoint& b) { + return Landpoint{ + a.lat * b.ati - a.ati * b.lat,//x + a.ati * b.lon - a.lon * b.ati,//y + a.lon * b.lat - a.lat * b.lon//z + }; +} + +float LAMPTOOLAPI cross2d(Point_3d a, Point_3d b) +{ + return a.x * b.y - a.y * b.x; +} + +Point_3d LAMPTOOLAPI operator -(Point_3d a, Point_3d b) +{ + return Point_3d{ a.x - b.x, a.y - b.y, a.z - b.z }; +} + +Point_3d LAMPTOOLAPI operator +(Point_3d a, Point_3d b) +{ + return Point_3d{ a.x + b.x, a.y + b.y, a.z + b.z }; +} + +double LAMPTOOLAPI operator /(Point_3d a, Point_3d b) +{ + return sqrt(pow(a.x, 2) + pow(a.y, 2)) / sqrt(pow(b.x, 2) + pow(b.y, 2)); +} + + +Landpoint LAMPTOOLAPI getSlopeVector(const Landpoint& p0, const Landpoint& p1, const Landpoint& p2, const Landpoint& p3, const Landpoint& p4) { + + Landpoint n0 = LLA2XYZ(p0), + n1 = LLA2XYZ(p1), + n2 = LLA2XYZ(p2), + n3 = LLA2XYZ(p3), + n4 = LLA2XYZ(p4); + Landpoint n01 = n1 - n0, n02 = n2 - n0, n03 = n3 - n0, n04 = n4 - n0; + // 锟斤拷n01为锟斤拷锟斤拷锟斤拷锟斤拷 + Landpoint np01 = p1 - p0, np02 = p2 - p0, np03 = p3 - p0, np04 = p4 - p0; + double a2 = getAngle(Landpoint{ np01.lon,np01.lat,0 }, Landpoint{ np02.lon,np02.lat,0 });// 01->02 锟斤拷时锟斤拷 + double a3 = getAngle(Landpoint{ np01.lon,np01.lat,0 }, Landpoint{ np03.lon,np03.lat,0 });// 01->03 锟斤拷时锟斤拷 + double a4 = getAngle(Landpoint{ np01.lon,np01.lat,0 }, Landpoint{ np04.lon,np04.lat,0 });// 01->04 锟斤拷时锟斤拷 + //qDebug() << a2 << "\t" << a3 << "\t" << a4 << endl; + a2 = 360 - a2; + a3 = 360 - a3; + a4 = 360 - a4; + Landpoint N, W, S, E; + N = n01; + if (a2 >= a3 && a2 >= a4) { + W = n02; + if (a3 >= a4) { + S = n03; + E = n04; + } + else { + S = n04; + E = n03; + } + } + else if (a3 >= a2 && a3 >= a4) { + W = n03; + if (a2 >= a4) { + S = n02; + E = n04; + } + else { + S = n04; + E = n02; + } + } + else if (a4 >= a2 && a4 >= a3) + { + W = n04; + if (a2 >= a3) { + S = n02; + E = n03; + } + else { + S = n03; + E = n02; + } + } + return (crossProduct(N, W) + crossProduct(W, S) + crossProduct(S, E) + crossProduct(E, N)) * 0.25; +} + + + +double LAMPTOOLAPI distance(const Vector3D& p1, const Vector3D& p2) +{ + double dx = p1.x - p2.x; + double dy = p1.y - p2.y; + double dz = p1.z - p2.z; + return std::sqrt(dx * dx + dy * dy + dz * dz); +} + +double LAMPTOOLAPI pointToLineDistance(const Vector3D& point, const Vector3D& linePoint, const Vector3D& lineDirection) +{ + Vector3D pointToLine = { point.x - linePoint.x, point.y - linePoint.y, point.z - linePoint.z }; + + // 璁＄畻鐐瑰埌鐩寸嚎鐨勬姇褰辩偣鐨勪綅缃� + double t = (pointToLine.x * lineDirection.x + pointToLine.y * lineDirection.y + pointToLine.z * lineDirection.z) / + (lineDirection.x * lineDirection.x + lineDirection.y * lineDirection.y + lineDirection.z * lineDirection.z); + + // 璁＄畻鎶曞奖鐐� + Vector3D projection = { linePoint.x + t * lineDirection.x, linePoint.y + t * lineDirection.y, linePoint.z + t * lineDirection.z }; + + // 璁＄畻鐐瑰埌鐩寸嚎鐨勮窛绂� + return distance(point, projection); +} + +SphericalCoordinates LAMPTOOLAPI cartesianToSpherical(const CartesianCoordinates& cartesian) +{ + SphericalCoordinates spherical; + + spherical.r = std::sqrt(cartesian.x * cartesian.x + cartesian.y * cartesian.y + cartesian.z * cartesian.z); + spherical.theta = std::acos(cartesian.z / spherical.r); + spherical.phi = std::atan2(cartesian.y, cartesian.x); + + return spherical; +} + +CartesianCoordinates LAMPTOOLAPI sphericalToCartesian(const SphericalCoordinates& spherical) +{ + CartesianCoordinates cartesian; + + cartesian.x = spherical.r * std::sin(spherical.theta) * std::cos(spherical.phi); + cartesian.y = spherical.r * std::sin(spherical.theta) * std::sin(spherical.phi); + cartesian.z = spherical.r * std::cos(spherical.theta); + + return cartesian; +} \ No newline at end of file diff --git a/src/LAMPTool/GeoOperator.h b/src/LAMPTool/GeoOperator.h new file mode 100644 index 0000000..6353971 --- /dev/null +++ b/src/LAMPTool/GeoOperator.h @@ -0,0 +1,111 @@ +#pragma once + + +#ifndef _GEOOPERATOR_H +#define _GEOOPERATOR_H +#include "LAMPToolAPI.h" +#include "BaseConstVariable.h" +#include +#include +#include +#include +#include +#include +#include + +/// +/// 三维向量，坐标表达 +/// +struct Landpoint // 点 SAR影像的像素坐标； +{ + /// + /// 经度x + /// + double lon; // 经度x lon pixel_col + /// + /// 纬度y + /// + double lat; // 纬度y lat pixel_row + /// + /// 高度z + /// + double ati; // 高程z ati pixel_time +}; +struct Point_3d { + double x; + double y; + double z; +}; + +/// +/// 将经纬度转换为地固参心坐标系 +/// +/// 经纬度点--degree +/// 投影坐标系点 +Landpoint LAMPTOOLAPI LLA2XYZ(const Landpoint& LLA); +Eigen::MatrixXd LAMPTOOLAPI LLA2XYZ(Eigen::MatrixXd landpoint); + +/// +/// 将地固参心坐标系转换为经纬度 +/// +/// 固参心坐标系 +/// 经纬度--degree +Landpoint LAMPTOOLAPI XYZ2LLA(const Landpoint& XYZ); + + +Landpoint LAMPTOOLAPI operator +(const Landpoint& p1, const Landpoint& p2); + +Landpoint LAMPTOOLAPI operator -(const Landpoint& p1, const Landpoint& p2); + +bool LAMPTOOLAPI operator ==(const Landpoint& p1, const Landpoint& p2); + +Landpoint LAMPTOOLAPI operator *(const Landpoint& p, double scale); + +double LAMPTOOLAPI getAngle(const Landpoint& a, const Landpoint& b); + +double LAMPTOOLAPI dot(const Landpoint& p1, const Landpoint& p2); + +double LAMPTOOLAPI getlength(const Landpoint& p1); + +Landpoint LAMPTOOLAPI crossProduct(const Landpoint& a, const Landpoint& b); + + +Landpoint LAMPTOOLAPI getSlopeVector(const Landpoint& p0, const Landpoint& p1, const Landpoint& p2, const Landpoint& p3, const Landpoint& p4); + + + +float LAMPTOOLAPI cross2d(Point_3d a, Point_3d b); + +Point_3d LAMPTOOLAPI operator -(Point_3d a, Point_3d b); + +Point_3d LAMPTOOLAPI operator +(Point_3d a, Point_3d b); + +double LAMPTOOLAPI operator /(Point_3d a, Point_3d b); + + + +// 矢量计算 +struct Vector3D { + double x, y, z; +}; + +// 计算两点之间的距离 +double LAMPTOOLAPI distance(const Vector3D& p1, const Vector3D& p2); +// 计算点到直线的最短距离 +double LAMPTOOLAPI pointToLineDistance(const Vector3D& point, const Vector3D& linePoint, const Vector3D& lineDirection); + + +struct CartesianCoordinates { + double x, y, z; +}; + +struct SphericalCoordinates { + double r, theta, phi; +}; + +SphericalCoordinates LAMPTOOLAPI cartesianToSpherical(const CartesianCoordinates& cartesian); + +CartesianCoordinates LAMPTOOLAPI sphericalToCartesian(const SphericalCoordinates& spherical); + + +#endif \ No newline at end of file diff --git a/src/LAMPTool/ImageOperatorBase.cpp b/src/LAMPTool/ImageOperatorBase.cpp new file mode 100644 index 0000000..9e75883 --- /dev/null +++ b/src/LAMPTool/ImageOperatorBase.cpp @@ -0,0 +1,1237 @@ +#include "ImageOperatorBase.h" +#include "GeoOperator.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "FileOperator.h" + +using namespace std; +using namespace Eigen; + +/** +* 输入数据是ENVI格式数据 +*/ + + + + + + +std::shared_ptr LAMPTOOLAPI OpenDataset(const QString& in_path,GDALAccess rwmode) +{ + GDALAllRegister(); + GDALDataset* dataset_ptr = (GDALDataset*)(GDALOpen(in_path.toUtf8().constData(), rwmode)); + std::shared_ptr rasterDataset(dataset_ptr, CloseDataset); + return rasterDataset; +} + +void LAMPTOOLAPI CloseDataset(GDALDataset* ptr) +{ + GDALClose(ptr); + ptr = NULL; +} + +int LAMPTOOLAPI TIFF2ENVI(QString in_tiff_path, QString out_envi_path) +{ + std::shared_ptr ds = OpenDataset(in_tiff_path); + const char* args[] = { "-of", "ENVI", NULL }; + GDALTranslateOptions* psOptions = GDALTranslateOptionsNew((char**)args, NULL); + GDALClose(GDALTranslate(out_envi_path.toUtf8().constData(),ds.get(), psOptions,NULL)); + GDALTranslateOptionsFree(psOptions); + return 0; +} + +int LAMPTOOLAPI ENVI2TIFF(QString in_envi_path, QString out_tiff_path) +{ + std::shared_ptr ds = OpenDataset(in_envi_path); + const char* args[] = { "-of", "Gtiff", NULL }; + GDALTranslateOptions* psOptions = GDALTranslateOptionsNew((char**)args, NULL); + GDALClose(GDALTranslate(out_tiff_path.toUtf8().constData(), ds.get(), psOptions, NULL)); + GDALTranslateOptionsFree(psOptions); + return 0; +} + +int LAMPTOOLAPI CreateDataset(QString new_file_path,int height, int width, int band_num,double* gt, QString projection, GDALDataType gdal_dtype ,bool need_gt) +{ + GDALAllRegister(); + GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI"); + std::shared_ptr< GDALDataset> poDstDS(poDriver->Create(new_file_path.toUtf8().constData(), width, height, band_num, gdal_dtype, NULL)); + if (need_gt) + { + poDstDS->SetProjection(projection.toUtf8().constData()); + poDstDS->SetGeoTransform(gt); + } + else {} + GDALFlushCache((GDALDatasetH)poDstDS.get()); + return 0; +} + +int LAMPTOOLAPI saveDataset(QString new_file_path, int start_line,int start_cols ,int band_ids, int datacols,int datarows,void* databuffer) +{ + GDALAllRegister(); + std::shared_ptr poDstDS = OpenDataset(new_file_path,GA_Update); + GDALDataType gdal_datatype = poDstDS->GetRasterBand(1)->GetRasterDataType(); + poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_cols, start_line, datacols, datarows, databuffer, datacols, datarows, gdal_datatype, 0, 0); + GDALFlushCache(poDstDS.get()); + return 0; +} + +int LAMPTOOLAPI block_num_pre_memory(int block_width, int height, GDALDataType gdal_datatype,double memey_size) +{ + // 计算大小 + int size_meta = 0; + + if (gdal_datatype == GDT_Byte) { + size_meta = 1; + } + else if (gdal_datatype == GDT_UInt16) { + size_meta = 2; // 只有双通道才能构建 复数矩阵 + } + else if (gdal_datatype == GDT_UInt16) { + size_meta = 2; + } + else if (gdal_datatype == GDT_Int16) { + size_meta = 2; + } + else if (gdal_datatype == GDT_UInt32) { + size_meta = 4; + } + else if (gdal_datatype == GDT_Int32) { + size_meta = 4; + } + //else if (gdal_datatype == GDT_UInt64) { + // size_meta = 8; + //} + //else if (gdal_datatype == GDT_Int64) { + // size_meta = 8; + //} + else if (gdal_datatype == GDT_Float32) { + size_meta = 4; + } + else if (gdal_datatype == GDT_Float64) { + size_meta = 4; + } + else if (gdal_datatype == GDT_CInt16) { size_meta = 2; } + else if (gdal_datatype == GDT_CInt32) { size_meta = 2; } + else if (gdal_datatype == GDT_CFloat32) { size_meta = 4; } + else if (gdal_datatype == GDT_CFloat64) { size_meta = 8; } + else {} + int block_num = int(memey_size / (size_meta * block_width)); + block_num = block_num > height ? height : block_num; // 行数 + block_num = block_num < 1 ? 1 : block_num; + return block_num; +} + +Eigen::Matrix LAMPTOOLAPI ReadComplexMatrixData(int start_line, int width, int line_num, std::shared_ptr rasterDataset, GDALDataType gdal_datatype) +{ + + int band_num = rasterDataset->GetRasterCount(); + if (gdal_datatype == 0) { + return Eigen::Matrix (0, 0); + } + else if (gdal_datatype < 8) { + if (band_num != 2) { + return Eigen::Matrix(0, 0); + } + } + else if (gdal_datatype < 12) { + if (band_num != 1) { + return Eigen::Matrix(0, 0); + } + + } + else {} + bool _flag = false; + Eigen::Matrix data_mat(line_num * width, 2);// 必须强制行优先 + if (gdal_datatype == GDT_Byte) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX imag_mat(line_num * width, 1); + rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag + data_mat.col(0) = (real_mat.array().cast()).array(); + data_mat.col(1) = (imag_mat.array().cast()).array(); + _flag = true; + } + else if (gdal_datatype == GDT_UInt16) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX imag_mat(line_num * width, 1); + rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag + data_mat.col(0) = (real_mat.array().cast()).array(); + data_mat.col(1) = (imag_mat.array().cast()).array(); + _flag = true; + } + else if (gdal_datatype == GDT_Int16) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX imag_mat(line_num * width, 1); + rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag + data_mat.col(0) = (real_mat.array().cast()).array(); + data_mat.col(1) = (imag_mat.array().cast()).array(); + _flag = true; + } + else if (gdal_datatype == GDT_UInt32) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX imag_mat(line_num * width, 1); + rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag + data_mat.col(0) = (real_mat.array().cast()).array(); + data_mat.col(1) = (imag_mat.array().cast()).array(); + _flag = true; + } + else if (gdal_datatype == GDT_Int32) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX imag_mat(line_num * width, 1); + rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag + data_mat.col(0) = (real_mat.array().cast()).array(); + data_mat.col(1) = (imag_mat.array().cast()).array(); + _flag = true; + } + //else if (gdal_datatype == GDT_UInt64) { + // Eigen::MatrixX real_mat(line_num * width, 1); + // Eigen::MatrixX imag_mat(line_num * width, 1); + // rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + // rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag + // data_mat.col(0) = (real_mat.array().cast()).array(); + // data_mat.col(1) = (imag_mat.array().cast()).array(); + // _flag = true; + //} + //else if (gdal_datatype == GDT_Int64) { + // Eigen::MatrixX real_mat(line_num * width, 1); + // Eigen::MatrixX imag_mat(line_num * width, 1); + // rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + // rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag + // data_mat.col(0) = (real_mat.array().cast()).array(); + // data_mat.col(1) = (imag_mat.array().cast()).array(); + // _flag = true; + //} + else if (gdal_datatype == GDT_Float32) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX imag_mat(line_num * width, 1); + rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag + data_mat.col(0) = (real_mat.array().cast()).array(); + data_mat.col(1) = (imag_mat.array().cast()).array(); + _flag = true; + } + else if (gdal_datatype == GDT_Float64) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX imag_mat(line_num * width, 1); + rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num, imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag + data_mat.col(0) = (real_mat.array().cast()).array(); + data_mat.col(1) = (imag_mat.array().cast()).array(); + _flag = true; + } + else if (gdal_datatype == GDT_CInt16) { + Eigen::MatrixX> complex_short_mat(line_num * width, 1); + rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, 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()).array(); + data_mat.col(1) = (complex_short_mat.imag().array().cast()).array(); + _flag = true; + } + else if (gdal_datatype == GDT_CInt32) { + Eigen::MatrixX> complex_short_mat(line_num * width, 1); + rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, 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()).array(); + data_mat.col(1) = (complex_short_mat.imag().array().cast()).array(); + _flag = true; + } + else if (gdal_datatype == GDT_CFloat32) { + Eigen::MatrixX> complex_short_mat(line_num * width, 1); + rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, 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()).array(); + data_mat.col(1) = (complex_short_mat.imag().array().cast()).array(); + _flag = true; + } + else if (gdal_datatype == GDT_CFloat64) { + Eigen::MatrixX> complex_short_mat(line_num * width, 1); + rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, 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()).array(); + data_mat.col(1) = (complex_short_mat.imag().array().cast()).array(); + _flag = true; + } + else {} + // 保存数据 + + + if (_flag) { + return data_mat; + } + else { + return Eigen::Matrix(0, 0);// 必须强制行优先; + } +} + +Eigen::Matrix LAMPTOOLAPI ReadMatrixDoubleData(int start_line, int width, int line_num, std::shared_ptr rasterDataset, GDALDataType gdal_datatype, int band_idx) +{ + // 构建矩阵块，使用eigen 进行矩阵计算，加速计算 + bool _flag = false; + Eigen::Matrix data_mat(line_num * width, 1);// 必须强制行优先 + if (gdal_datatype == GDT_Byte) { + Eigen::MatrixX real_mat(line_num * width, 1); + rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; + _flag = true; + } + else if (gdal_datatype == GDT_UInt16) { + Eigen::MatrixX real_mat(line_num * width, 1); + rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; + _flag = true; + } + else if (gdal_datatype == GDT_Int16) { + Eigen::MatrixX real_mat(line_num * width, 1); + rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; + _flag = true; + } + else if (gdal_datatype == GDT_UInt32) { + Eigen::MatrixX real_mat(line_num * width, 1); + rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; + _flag = true; + } + else if (gdal_datatype == GDT_Int32) { + Eigen::MatrixX real_mat(line_num * width, 1); + rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; + _flag = true; + } + //else if (gdal_datatype == GDT_UInt64) { + // Eigen::MatrixX real_mat(line_num * width, 1); + // rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + // data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; + // _flag = true; + //} + //else if (gdal_datatype == GDT_Int64) { + // Eigen::MatrixX real_mat(line_num * width, 1); + // rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + // data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; + // _flag = true; + //} + else if (gdal_datatype == GDT_Float32) { + Eigen::MatrixX real_mat(line_num * width, 1); + rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; + _flag = true; + } + else if (gdal_datatype == GDT_Float64) { + Eigen::MatrixX real_mat(line_num * width, 1); + rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num, real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real + data_mat.col(0) = ((real_mat.array().cast()).array().pow(2)).log10() * 10.0; + _flag = true; + } + else {} + + + return data_mat; +} + +Eigen::MatrixXd LAMPTOOLAPI getGeoTranslationArray(QString in_path) +{ + + + + return Eigen::MatrixXd(); +} + +ImageGEOINFO LAMPTOOLAPI getImageINFO(QString in_path) +{ + std::shared_ptr df = OpenDataset(in_path); + int width = df->GetRasterXSize(); + int heigh = df->GetRasterYSize(); + int band_num = df->GetRasterCount(); + ImageGEOINFO result; + result.width = width; + result.height = heigh; + result.bandnum = band_num; + + return result; +} + +GDALDataType LAMPTOOLAPI getGDALDataType(QString fileptah) +{ + omp_lock_t lock; + omp_init_lock(&lock); + omp_set_lock(&lock); + GDALAllRegister(); + + GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(fileptah.toUtf8().constData(), GA_ReadOnly));//锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷 + + GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType(); + + GDALClose((GDALDatasetH)rasterDataset); + omp_unset_lock(&lock); //锟酵放伙拷斤拷 + omp_destroy_lock(&lock); //劫伙拷斤拷 + + return gdal_datatype; +} + + +gdalImage::gdalImage() +{ + this->height = 0; + this->width = 0; + this->data_band_ids = 1; + this->start_row = 0; + this->start_col = 0; +} + +/// +/// 斤拷图饺∮帮拷锟?1锟?7 +/// +/// +gdalImage::gdalImage(const QString& raster_path) +{ + omp_lock_t lock; + omp_init_lock(&lock); // 锟斤拷始斤拷斤拷 + omp_set_lock(&lock); //锟斤拷没斤拷锟?1锟?7 + this->img_path = raster_path; + + GDALAllRegister();// 注绞斤拷斤拷锟?1锟?7 + // 锟斤拷DEM影锟斤拷 + GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(raster_path.toUtf8().constData(), GA_ReadOnly));//锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷 + this->width = rasterDataset->GetRasterXSize(); + this->height = rasterDataset->GetRasterYSize(); + this->band_num = rasterDataset->GetRasterCount(); + + double* gt = new double[6]; + // 锟斤拷梅斤拷锟斤拷 + rasterDataset->GetGeoTransform(gt); + this->gt = Eigen::MatrixXd(2, 3); + this->gt << gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]; + + this->projection = rasterDataset->GetProjectionRef(); + // 斤拷斤拷 + //double* inv_gt = new double[6];; + //GDALInvGeoTransform(gt, inv_gt); // 斤拷斤拷 + // 斤拷投影 + GDALFlushCache((GDALDatasetH)rasterDataset); + GDALClose((GDALDatasetH)rasterDataset); + rasterDataset = NULL;// 指矫匡拷 + this->InitInv_gt(); + delete[] gt; + ////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH + omp_unset_lock(&lock); //锟酵放伙拷斤拷 + omp_destroy_lock(&lock); //劫伙拷斤拷 +} + +gdalImage::~gdalImage() +{ +} + +void gdalImage::setHeight(int height) +{ + this->height = height; +} + +void gdalImage::setWidth(int width) +{ + this->width = width; +} + +void gdalImage::setTranslationMatrix(Eigen::MatrixXd gt) +{ + this->gt = gt; +} + +void gdalImage::setData(Eigen::MatrixXd, int data_band_ids) +{ + this->data = data; + this->data_band_ids = data_band_ids; +} + + + +Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count, int cols_count, int band_ids = 1) +{ + omp_lock_t lock; + omp_init_lock(&lock); + omp_set_lock(&lock); + GDALAllRegister(); + + GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly));//锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷 + + GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType(); + GDALRasterBand* demBand = rasterDataset->GetRasterBand(band_ids); + + rows_count = start_row + rows_count <= this->height ? rows_count : this->height - start_row; + cols_count = start_col + cols_count <= this->width ? cols_count : this->width - start_col; + + MatrixXd datamatrix(rows_count, cols_count); + + + if (gdal_datatype == GDT_Byte) { + char* temp = new char[rows_count * cols_count]; + 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]; + } + } + delete[] temp; + } + else if (gdal_datatype == GDT_UInt16) { + unsigned short* temp = new unsigned short[rows_count * cols_count]; + 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]; + } + } + delete[] temp; + } + else if (gdal_datatype == GDT_Int16) { + short* temp = new short[rows_count * cols_count]; + 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]; + } + } + delete[] temp; + } + else if (gdal_datatype == GDT_UInt32) { + unsigned int* temp = new unsigned int[rows_count * cols_count]; + 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]; + } + } + delete[] temp; + } + else if (gdal_datatype == GDT_Int32) { + int* temp = new int[rows_count * cols_count]; + 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]; + } + } + delete[] temp; + } + //else if (gdal_datatype == GDT_UInt64) { + // unsigned long* temp = new unsigned long[rows_count * cols_count]; + // 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]; + // } + // } + // delete[] temp; + //} + //else if (gdal_datatype == GDT_Int64) { + // long* temp = new long[rows_count * cols_count]; + // 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]; + // } + // } + // delete[] temp; + //} + else if (gdal_datatype == GDT_Float32) { + float* temp = new float[rows_count * cols_count]; + 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]; + } + } + delete[] temp; + } + else if (gdal_datatype == GDT_Float64) { + double* temp = new double[rows_count * cols_count]; + 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]; + } + } + delete[] temp; + } + else {} + GDALClose((GDALDatasetH)rasterDataset); + omp_unset_lock(&lock); //锟酵放伙拷斤拷 + omp_destroy_lock(&lock); //劫伙拷斤拷 + //GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH + return datamatrix; + +} + +Eigen::MatrixXd gdalImage::getGeoTranslation() +{ + return this->gt; +} + +GDALDataType gdalImage::getDataType() +{ + GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly)); + GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType(); + return gdal_datatype; +} + +/// +/// 写锟斤拷遥锟斤拷影锟斤拷 +/// +/// +/// +/// +/// +void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col = 0, int band_ids = 1) +{ + omp_lock_t lock; + omp_init_lock(&lock); + omp_set_lock(&lock); + if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) { + QString tip = u8"file path: " + this->img_path; + qDebug() << tip; + throw exception(tip.toUtf8().constData()); + } + GDALAllRegister(); + GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); + GDALDataset* poDstDS = nullptr; + if (exists_test(this->img_path)) { + poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update)); + } + else { + poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height, this->band_num, GDT_Float32, NULL); // 斤拷锟斤拷 + poDstDS->SetProjection(this->projection.toUtf8().constData()); + + + double gt_ptr[6]; + for (int i = 0; i < this->gt.rows(); i++) { + for (int j = 0; j < this->gt.cols(); j++) { + gt_ptr[i * 3 + j] = this->gt(i, j); + } + } + poDstDS->SetGeoTransform(gt_ptr); + delete[] gt_ptr; + } + + 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; + } + } + //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, GDT_Float32, 0, 0); + + GDALFlushCache(poDstDS); + GDALClose((GDALDatasetH)poDstDS); + //GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH + delete[] databuffer; + omp_unset_lock(&lock); //锟酵放伙拷斤拷 + omp_destroy_lock(&lock); //劫伙拷斤拷 +} + +void gdalImage::saveImage() +{ + this->saveImage(this->data, this->start_row, this->start_col, this->data_band_ids); +} + +void gdalImage::setNoDataValue(double nodatavalue = -9999, int band_ids = 1) +{ + GDALAllRegister();// 注绞斤拷斤拷锟?1锟?7 + //GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); + GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.toUtf8().constData(), GA_Update)); + poDstDS->GetRasterBand(band_ids)->SetNoDataValue(nodatavalue); + GDALFlushCache((GDALDatasetH)poDstDS); + GDALClose((GDALDatasetH)poDstDS); +} + +int gdalImage::InitInv_gt() +{ + //1 lon lat = x + //1 lon lat = y + Eigen::MatrixXd temp=Eigen::MatrixXd::Zero(2, 3); + this->inv_gt = temp; + double a = this->gt(0, 0); + double b = this->gt(0, 1); + double c = this->gt(0, 2); + double d = this->gt(1, 0); + double e = this->gt(1, 1); + double f = this->gt(1, 2); + double g = 1; + double det_gt = b * f - c * e; + if (det_gt == 0) { + return 0; + } + this->inv_gt(0, 0) = (c * d - a * f) / det_gt; //2 + this->inv_gt(0, 1) = f / det_gt; // lon + this->inv_gt(0, 2) = -c / det_gt; // lat + this->inv_gt(1, 0) = (a * e - b * d) / det_gt; //1 + this->inv_gt(1, 1) = -e / det_gt; // lon + this->inv_gt(1, 2) = b / det_gt; // lat + return 1; +} + +Landpoint gdalImage::getRow_Col(double lon, double lat) +{ + Landpoint p{ 0,0,0 }; + p.lon = this->inv_gt(0, 0) + lon * this->inv_gt(0, 1) + lat * this->inv_gt(0, 2); //x + p.lat = this->inv_gt(1, 0) + lon * this->inv_gt(1, 1) + lat * this->inv_gt(1, 2); //y + return p; +} + +Landpoint gdalImage::getLandPoint(double row, double col, double ati = 0) +{ + Landpoint p{ 0,0,0 }; + p.lon = this->gt(0, 0) + col * this->gt(0, 1) + row * this->gt(0, 2); //x + p.lat = this->gt(1, 0) + col * this->gt(1, 1) + row * this->gt(1, 2); //y + p.ati = ati; + return p; +} + +double gdalImage::mean(int bandids) +{ + double mean_value = 0; + double count = this->height * this->width; + int line_invert = 100; + int start_ids = 0; + do { + Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids); + mean_value = mean_value + sar_a.sum() / count; + start_ids = start_ids + line_invert; + } while (start_ids < this->height); + return mean_value; +} + +double gdalImage::max(int bandids) +{ + double max_value = 0; + bool state_max = true; + int line_invert = 100; + int start_ids = 0; + double temp_max = 0; + do { + Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids); + if (state_max) { + state_max = false; + max_value = sar_a.maxCoeff(); + } + else { + temp_max = sar_a.maxCoeff(); + if (max_value < temp_max) { + max_value = temp_max; + } + } + start_ids = start_ids + line_invert; + } while (start_ids < this->height); + return max_value; +} + +double gdalImage::min(int bandids) +{ + double min_value = 0; + bool state_min = true; + int line_invert = 100; + int start_ids = 0; + double temp_min = 0; + do { + Eigen::MatrixXd sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids); + if (state_min) { + state_min = false; + min_value = sar_a.minCoeff(); + } + else { + temp_min = sar_a.minCoeff(); + if (min_value < temp_min) { + min_value = temp_min; + } + } + start_ids = start_ids + line_invert; + } while (start_ids < this->height); + return min_value; +} + +GDALRPCInfo gdalImage::getRPC() +{ + CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); + CPLSetConfigOption("GDAL_DATA", "./data"); + GDALAllRegister();//注斤拷锟斤拷 + //斤拷锟斤拷 + GDALDataset* pDS = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly); + //锟斤拷元斤拷锟叫伙拷取RPC锟斤拷息 + char** papszRPC = pDS->GetMetadata("RPC"); + + //斤拷取锟斤拷RPC锟斤拷息斤拷山峁癸拷锟?1锟?7 + GDALRPCInfo oInfo; + GDALExtractRPCInfo(papszRPC, &oInfo); + + GDALClose((GDALDatasetH)pDS); + + return oInfo; +} + +Eigen::MatrixXd gdalImage::getLandPoint(Eigen::MatrixXd points) +{ + if (points.cols() != 3) { + throw new exception("the size of points is equit 3!!!"); + } + + Eigen::MatrixXd result(points.rows(), 3); + result.col(2) = points.col(2);// 锟竭筹拷 + points.col(2) = points.col(2).array() * 0 + 1; + points.col(0).swap(points.col(2));// 斤拷 + Eigen::MatrixXd gts(3, 2); + gts.col(0) = this->gt.row(0); + gts.col(1) = this->gt.row(1); + + result.block(0, 0, points.rows(), 2) = points * gts; + return result; +} + +Eigen::MatrixXd gdalImage::getHist(int bandids) +{ + GDALAllRegister();// 注绞斤拷斤拷锟?1锟?7 + // 锟斤拷DEM影锟斤拷 + GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly));//锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷 + + GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType(); + GDALRasterBand* xBand = rasterDataset->GetRasterBand(bandids); + + + double dfMin = this->min(bandids); + double dfMax = this->max(bandids); + int count = int((dfMax - dfMin) / 0.01); + count = count > 255 ? count : 255; + GUIntBig* panHistogram = new GUIntBig[count]; + xBand->GetHistogram(dfMin, dfMax, count, panHistogram, TRUE, FALSE, NULL, NULL); + Eigen::MatrixXd result(count, 2); + double delta = (dfMax - dfMin) / count; + for (int i = 0; i < count; i++) { + result(i, 0) = dfMin + i * delta; + result(i, 1) = double(panHistogram[i]); + } + delete[] panHistogram; + GDALClose((GDALDatasetH)rasterDataset); + return result; +} + + +gdalImage LAMPTOOLAPI CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt,bool overwrite) { + + 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();// 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7 + GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); + GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, GDT_Float32, NULL); // 锟斤拷锟斤拷锟斤拷 + if (need_gt) { + 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; +} + +gdalImageComplex LAMPTOOLAPI CreategdalImageComplex(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite) +{ + if (exists_test(img_path.toUtf8().constData())) { + if (overwrite) { + gdalImageComplex result_img(img_path); + return result_img; + } + else { + throw "file has exist!!!"; + exit(1); + } + } + GDALAllRegister(); + GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI"); + GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, GDT_CFloat64, NULL); + if (need_gt) { + 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 + gdalImageComplex result_img(img_path); + return result_img; +} + + +int LAMPTOOLAPI ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int new_width, int new_height, GDALResampleAlg eResample) +{ + GDALAllRegister(); + CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); + GDALDataset* pDSrc = (GDALDataset*)GDALOpen(pszSrcFile, GA_ReadOnly); + if (pDSrc == NULL) + { + return -1; + } + + GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("GTiff"); + if (pDriver == NULL) + { + GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc); + return -2; + } + int width = pDSrc->GetRasterXSize(); + int height = pDSrc->GetRasterYSize(); + int nBandCount = pDSrc->GetRasterCount(); + GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType(); + + char* pszSrcWKT = NULL; + pszSrcWKT = const_cast(pDSrc->GetProjectionRef()); + + //锟斤拷锟矫伙拷锟酵队帮拷锟斤拷锟轿?拷锟斤拷锟揭伙拷锟?1锟?7 + if (strlen(pszSrcWKT) <= 0) + { + OGRSpatialReference oSRS; + oSRS.importFromEPSG(4326); + //oSRS.SetUTM(50, true); //锟斤拷锟斤拷锟斤拷 锟斤拷锟斤拷120锟斤拷 + //oSRS.SetWellKnownGeogCS("WGS84"); + oSRS.exportToWkt(&pszSrcWKT); + } + qDebug() << "GDALCreateGenImgProjTransformer " << endl; + void* hTransformArg; + hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT, FALSE, 0.0, 1); + qDebug() << "no proj " << endl; + //(没锟斤拷投影锟斤拷影锟斤拷锟斤拷锟斤拷锟斤拷卟锟酵?拷锟?1锟?7) + if (hTransformArg == NULL) + { + GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc); + return -3; + } + qDebug() << "has proj " << endl; + double dGeoTrans[6] = { 0 }; + int nNewWidth = 0, nNewHeight = 0; + if (GDALSuggestedWarpOutput((GDALDatasetH)pDSrc, GDALGenImgProjTransform, hTransformArg, dGeoTrans, &nNewWidth, &nNewHeight) != CE_None) + { + GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc); + return -3; + } + + //GDALDestroyGenImgProjTransformer(hTransformArg); + + + GDALDataset* pDDst = pDriver->Create(pszOutFile, new_width, new_height, nBandCount, dataType, NULL); + if (pDDst == NULL) + { + GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc); + return -2; + } + + pDDst->SetProjection(pszSrcWKT); + pDDst->SetGeoTransform(gt); + + GDALWarpOptions* psWo = GDALCreateWarpOptions(); + + //psWo->papszWarpOptions = CSLDuplicate(NULL); + psWo->eWorkingDataType = dataType; + psWo->eResampleAlg = eResample; + + psWo->hSrcDS = (GDALDatasetH)pDSrc; + psWo->hDstDS = (GDALDatasetH)pDDst; + qDebug() << "GDALCreateGenImgProjTransformer" << endl; + psWo->pfnTransformer = GDALGenImgProjTransform; + psWo->pTransformerArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, (GDALDatasetH)pDDst, pszSrcWKT, FALSE, 0.0, 1);; + + qDebug() << "GDALCreateGenImgProjTransformer has created" << endl; + psWo->nBandCount = nBandCount; + psWo->panSrcBands = (int*)CPLMalloc(nBandCount * sizeof(int)); + psWo->panDstBands = (int*)CPLMalloc(nBandCount * sizeof(int)); + for (int i = 0; i < nBandCount; i++) + { + psWo->panSrcBands[i] = i + 1; + psWo->panDstBands[i] = i + 1; + } + + GDALWarpOperation oWo; + if (oWo.Initialize(psWo) != CE_None) + { + GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc); + GDALClose((GDALDatasetH)(GDALDatasetH)pDDst); + return -3; + } + qDebug() << "ChunkAndWarpImage:" << new_width << "," << new_height << endl; + oWo.ChunkAndWarpMulti(0, 0, new_width, new_height); + GDALFlushCache(pDDst); + qDebug() << "ChunkAndWarpImage over" << endl; + //GDALDestroyGenImgProjTransformer(psWo->pTransformerArg); + //GDALDestroyWarpOptions(psWo); + GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc); + GDALClose((GDALDatasetH)(GDALDatasetH)pDDst); + ////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH + return 0; +} + +int LAMPTOOLAPI ResampleGDALs(const char* pszSrcFile, int band_ids, GDALRIOResampleAlg eResample) +{ + GDALAllRegister(); + CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO"); + GDALDataset* pDSrc = (GDALDataset*)GDALOpen(pszSrcFile, GA_Update); + if (pDSrc == NULL) + { + return -1; + } + + GDALDataType gdal_datatype = pDSrc->GetRasterBand(1)->GetRasterDataType(); + + GDALRasterBand* demBand = pDSrc->GetRasterBand(band_ids); + + int width = pDSrc->GetRasterXSize(); + int height = pDSrc->GetRasterYSize(); + int start_col = 0, start_row = 0, rows_count = 0, cols_count; + + int row_delta = int(120000000 / width); + + GDALRasterIOExtraArg psExtraArg; + INIT_RASTERIO_EXTRA_ARG(psExtraArg); + psExtraArg.eResampleAlg = eResample; + + do { + + rows_count = start_row + row_delta < height ? row_delta : height - start_row; + cols_count = width; + + if (gdal_datatype == GDALDataType::GDT_UInt16) { + + unsigned short* temp = new unsigned short[rows_count * cols_count]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); + demBand->RasterIO(GF_Write, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0, &psExtraArg); + delete[] temp; + + } + else if (gdal_datatype == GDALDataType::GDT_Int16) { + + short* temp = new short[rows_count * cols_count]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); + demBand->RasterIO(GF_Write, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0, &psExtraArg); + delete[] temp; + } + else if (gdal_datatype == GDALDataType::GDT_Float32) { + float* temp = new float[rows_count * cols_count]; + demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0); + demBand->RasterIO(GF_Write, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0, &psExtraArg); + delete[] temp; + } + start_row = start_row + rows_count; + } while (start_row < height); + GDALClose((GDALDatasetH)pDSrc); + + + return 0; +} + +int LAMPTOOLAPI saveMatrixXcd2TiFF(Eigen::MatrixXcd data, QString out_tiff_path) +{ + int rows = data.rows(); + int cols = data.cols(); + + Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3); + + gdalImage image_tiff = CreategdalImage(out_tiff_path, rows, cols, 2, gt, "", false, true);// 注意这里保留仿真结果 + // 保存二进制文件 + Eigen::MatrixXd real_img = data.array().real(); + Eigen::MatrixXd imag_img = data.array().imag(); + image_tiff.saveImage(real_img, 0, 0, 1); + image_tiff.saveImage(imag_img, 0, 0, 2); + return -1; +} + + +gdalImageComplex::gdalImageComplex(const QString& raster_path) +{ + omp_lock_t lock; + omp_init_lock(&lock); + omp_set_lock(&lock); + this->img_path = raster_path; + + GDALAllRegister(); + GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(raster_path.toUtf8().constData(), GA_ReadOnly));//锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷 + this->width = rasterDataset->GetRasterXSize(); + this->height = rasterDataset->GetRasterYSize(); + this->band_num = rasterDataset->GetRasterCount(); + + double* gt = new double[6]; + rasterDataset->GetGeoTransform(gt); + this->gt = Eigen::MatrixXd(2, 3); + this->gt << gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]; + + + double a = this->gt(0, 0); + double b = this->gt(0, 1); + double c = this->gt(0, 2); + double d = this->gt(1, 0); + double e = this->gt(1, 1); + double f = this->gt(1, 2); + + this->projection = rasterDataset->GetProjectionRef(); + + // 斤拷投影 + GDALFlushCache((GDALDatasetH)rasterDataset); + GDALClose((GDALDatasetH)rasterDataset); + rasterDataset = NULL;// 指矫匡拷 + this->InitInv_gt(); + delete[] gt; + ////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH + omp_unset_lock(&lock); //锟酵放伙拷斤拷 + omp_destroy_lock(&lock); //劫伙拷斤拷 + + +} + +gdalImageComplex::~gdalImageComplex() +{ +} + +void gdalImageComplex::setData(Eigen::MatrixXcd data) +{ + this->data = data; +} + +void gdalImageComplex::saveImage(Eigen::MatrixXcd data, int start_row, int start_col, int band_ids) +{ + omp_lock_t lock; + omp_init_lock(&lock); + omp_set_lock(&lock); + if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) { + QString tip = u8"file path: " + this->img_path; + qDebug() << tip; + throw exception(tip.toUtf8().constData()); + } + GDALAllRegister(); + GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI"); + GDALDataset* poDstDS = nullptr; + if (exists_test(this->img_path)) { + poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update)); + } + else { + poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height, this->band_num, GDT_CFloat64, NULL); // 斤拷锟斤拷 + poDstDS->SetProjection(this->projection.toUtf8().constData()); + + + double gt_ptr[6]; + for (int i = 0; i < this->gt.rows(); i++) { + for (int j = 0; j < this->gt.cols(); j++) { + gt_ptr[i * 3 + j] = this->gt(i, j); + } + } + poDstDS->SetGeoTransform(gt_ptr); + delete[] gt_ptr; + } + + int datarows = data.rows(); + int datacols = data.cols(); + + double* databuffer = new double[data.size() * 2]; + for (int i = 0; i < data.rows(); i++) { + for (int j = 0; j < data.cols(); j++) { + databuffer[i * data.cols() * 2 + j * 2] = data(i, j).real(); + databuffer[i * data.cols() * 2 + j * 2 + 1] = data(i, j).imag(); + } + } + + //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, GDT_CFloat64, 0, 0); + GDALFlushCache(poDstDS); + delete databuffer; + GDALClose((GDALDatasetH)poDstDS); + //GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH + omp_unset_lock(&lock); //锟酵放伙拷斤拷 + omp_destroy_lock(&lock); //劫伙拷斤拷 + +} + +Eigen::MatrixXcd gdalImageComplex::getDataComplex(int start_row, int start_col, int rows_count, int cols_count, int band_ids) +{ + + + GDALDataset* poDataset; + GDALAllRegister(); + + // 打开TIFF文件 + poDataset = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly); + if (poDataset == nullptr) { + QMessageBox::warning(nullptr, u8"错误", u8"无法打开文件：" + this->img_path); + qDebug() << u8"无法打开文件：" + this->img_path; + } + + // 获取数据集的第一个波段 + GDALRasterBand* poBand; + poBand = poDataset->GetRasterBand(1); + + // 读取波段信息，假设是复数类型 + int nXSize = poBand->GetXSize(); + int nYSize = poBand->GetYSize(); + + + double* databuffer = new double[nXSize * nYSize * 2]; + poBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, databuffer, cols_count, rows_count, GDT_CFloat64, 0, 0); + GDALClose((GDALDatasetH)poDataset); + Eigen::MatrixXcd rasterData(nYSize, nXSize); // 使用Eigen的MatrixXcd + for (size_t i = 0; i < nYSize; i++) + { + for (size_t j = 0; j < nXSize; j++) { + rasterData(i, j) = std::complex(databuffer[i * nXSize * 2 + j * 2], databuffer[i * nXSize * 2 + j * 2 + 1]); + } + } + + delete databuffer; + return rasterData; + +} + +void gdalImageComplex::saveImage() +{ + this->saveImage(this->data, this->start_row, this->start_col, this->data_band_ids); +} \ No newline at end of file diff --git a/src/LAMPTool/ImageOperatorBase.h b/src/LAMPTool/ImageOperatorBase.h new file mode 100644 index 0000000..d731e38 --- /dev/null +++ b/src/LAMPTool/ImageOperatorBase.h @@ -0,0 +1,190 @@ +#pragma once +/** + * 影像数据操作项 + * 所有数据相关的操作，都是基于ENVI的数据格式 + * 影像读写操作时基于 GDAL 库 + * **/ + + + +#ifndef IMAGEOPERATORBASE_H +#define IMAGEOPERATORBASE_H + +#include "BaseConstVariable.h" +#include "GeoOperator.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // for CPLMalloc() +#include "referenceHeader.h" +using namespace std; +using namespace Eigen; + +struct ImageGEOINFO { + int width; + int height; + int bandnum; +}; + + +// 判断是否需要输出为DLL +#define DLLOUT +// 文件打开 +std::shared_ptr LAMPTOOLAPI OpenDataset(const QString& in_path, GDALAccess rwmode= GA_ReadOnly); // 当指令销毁时，调用GDALClose 销毁类型 +void LAMPTOOLAPI CloseDataset(GDALDataset* ptr); + +// 数据格式转换 +int LAMPTOOLAPI TIFF2ENVI(QString in_tiff_path,QString out_envi_path); +int LAMPTOOLAPI ENVI2TIFF(QString in_envi_path,QString out_tiff_path); + +// 保存影像数据 --直接保存 ENVI 文件 + +int LAMPTOOLAPI CreateDataset(QString new_file_path, int height, int width, int band_num, double* gt, QString projection, GDALDataType gdal_dtype, bool need_gt); // 创建文件 + +int LAMPTOOLAPI saveDataset(QString new_file_path, int start_line, int start_cols, int band_ids, int datacols, int datarows, void* databuffer); + +// 根据限制条件估算分块大小 +int LAMPTOOLAPI block_num_pre_memory(int width, int height, GDALDataType gdal_dtype,double memey_size); + +// 将结果转换为复数 或者 实数 +Eigen::Matrix LAMPTOOLAPI ReadComplexMatrixData(int start_line,int width, int line_num, std::shared_ptr rasterDataset, GDALDataType gdal_datatype); + +Eigen::Matrix LAMPTOOLAPI ReadMatrixDoubleData(int start_line, int width, int line_num, std::shared_ptr rasterDataset, GDALDataType gdal_datatype,int band_idx); + +Eigen::MatrixXd LAMPTOOLAPI getGeoTranslationArray(QString in_path); +ImageGEOINFO LAMPTOOLAPI getImageINFO(QString in_path); + +GDALDataType LAMPTOOLAPI getGDALDataType(QString fileptah); + + +struct DemBox { + double min_lat; //纬度 + double min_lon;//经度 + double max_lat;//纬度 + double max_lon;//经度 +}; + + +/// +/// gdalImage图像操作类 +/// +class LAMPTOOLAPI gdalImage +{ + +public: // 方法 + gdalImage(); + gdalImage(const QString& raster_path); + ~gdalImage(); + virtual void setHeight(int); + virtual void setWidth(int); + virtual void setTranslationMatrix(Eigen::MatrixXd gt); + virtual void setData(Eigen::MatrixXd,int data_band_ids=1); + virtual Eigen::MatrixXd getData(int start_row, int start_col, int rows_count, int cols_count, int band_ids); + virtual Eigen::MatrixXd getGeoTranslation(); + virtual GDALDataType getDataType(); + virtual void saveImage(Eigen::MatrixXd, int start_row, int start_col, int band_ids); + virtual void saveImage(); + virtual void setNoDataValue(double nodatavalue, int band_ids); + virtual int InitInv_gt(); + virtual Landpoint getRow_Col(double lon, double lat); + virtual Landpoint getLandPoint(double i, double j, double ati); + virtual double mean(int bandids = 1); + virtual double max(int bandids = 1); + virtual double min(int bandids = 1); + virtual GDALRPCInfo getRPC(); + virtual Eigen::MatrixXd getLandPoint(Eigen::MatrixXd points); + + virtual Eigen::MatrixXd getHist(int bandids); +public: + QString img_path; // 图像文件 + int height; // 高 + int width; // 宽 + int band_num;// 波段数 + int start_row;// + int start_col;// + int data_band_ids; + Eigen::MatrixXd gt; // 变换矩阵 + Eigen::MatrixXd inv_gt; // 逆变换矩阵 + Eigen::MatrixXd data; + QString projection; +}; + + + + + + + +/// +/// gdalImage图像操作类 +/// +class LAMPTOOLAPI gdalImageComplex:public gdalImage +{ + +public: // 方法 + gdalImageComplex(const QString& raster_path); + ~gdalImageComplex(); + void setData(Eigen::MatrixXcd); + void saveImage(Eigen::MatrixXcd data, int start_row, int start_col, int band_ids); + Eigen::MatrixXcd getDataComplex(int start_row, int start_col, int rows_count, int cols_count, int band_ids); + void saveImage() override; + +public: + Eigen::MatrixXcd data; +}; + + + + + + + + + + + + + + +gdalImage LAMPTOOLAPI CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false); +gdalImageComplex LAMPTOOLAPI CreategdalImageComplex(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt = true, bool overwrite = false); + + + + +int LAMPTOOLAPI ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int new_width, int new_height, GDALResampleAlg eResample); + +int LAMPTOOLAPI ResampleGDALs(const char* pszSrcFile, int band_ids, GDALRIOResampleAlg eResample = GRIORA_Bilinear); + + +//--------------------- 保存文博 ------------------------------- + +int LAMPTOOLAPI saveMatrixXcd2TiFF(Eigen::MatrixXcd data, QString out_tiff_path); + +//---------------------------------------------------- + + + +#ifndef DLLOUT + + + + + +#else +//#define DllExport __declspec( dllexport ) +//double __declspec(dllexport) ProcessMGCMathX_MGC(int Xbetaidx, int Xbwidx, double XTao, double satH, char* sigma_path, char* output_path, +// double p1_x, double p1_y, double p2_x, double p2_y, double p3_x, double p3_y, double p4_x, double p4_y) + +#endif + +#endif \ No newline at end of file diff --git a/src/LAMPTool/LAMPToolAPI.h b/src/LAMPTool/LAMPToolAPI.h new file mode 100644 index 0000000..f27183a --- /dev/null +++ b/src/LAMPTool/LAMPToolAPI.h @@ -0,0 +1,13 @@ +#ifndef _LAMPTOOLAPI_H_ +#define _LAMPTOOLAPI_H_ + +#include + + +#if defined(LAMPTOOL_API) +#define LAMPTOOLAPI Q_DECL_EXPORT +#else +#define LAMPTOOLAPI Q_DECL_IMPORT +#endif + +#endif diff --git a/src/LAMPTool/SARBaseTool.cpp b/src/LAMPTool/SARBaseTool.cpp new file mode 100644 index 0000000..d8b8e32 --- /dev/null +++ b/src/LAMPTool/SARBaseTool.cpp @@ -0,0 +1,283 @@ +// SARBaseTool.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。 +// + +#include "ImageOperatorBase.h" +#include "SARBaseTool.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + +int LAMPTOOLAPI Complex2dB(QString in_tiff, QString out_dB_path) +{ + GDALAllRegister(); + std::shared_ptr 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 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 data_mat(line_num * width, 1);// 必须强制行优先 + if (gdal_datatype == GDT_Byte) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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()).array().pow(2) + (imag_mat.array().cast()).array().pow(2)).log10() * 10.0; + _flag = true; + } + else if (gdal_datatype == GDT_UInt16) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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()).array().pow(2) + (imag_mat.array().cast()).array().pow(2)).log10() * 10.0; + _flag = true; + } + else if (gdal_datatype == GDT_Int16) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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()).array().pow(2) + (imag_mat.array().cast()).array().pow(2)).log10() * 10.0; + _flag = true; + } + else if (gdal_datatype == GDT_UInt32) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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()).array().pow(2) + (imag_mat.array().cast()).array().pow(2)).log10() * 10.0; + _flag = true; + } + else if (gdal_datatype == GDT_Int32) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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()).array().pow(2) + (imag_mat.array().cast()).array().pow(2)).log10() * 10.0; + _flag = true; + } + //else if (gdal_datatype == GDT_UInt64) { + // Eigen::MatrixX real_mat(line_num * width, 1); + // Eigen::MatrixX 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()).array().pow(2) + (imag_mat.array().cast()).array().pow(2)).log10() * 10.0; + // _flag = true; + //} + //else if (gdal_datatype == GDT_Int64) { + // Eigen::MatrixX real_mat(line_num * width, 1); + // Eigen::MatrixX 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()).array().pow(2) + (imag_mat.array().cast()).array().pow(2)).log10() * 10.0; + // _flag = true; + //} + else if (gdal_datatype == GDT_Float32) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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()).array().pow(2) + (imag_mat.array().cast()).array().pow(2)).log10() * 10.0; + _flag = true; + } + else if (gdal_datatype == GDT_Float64) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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()).array().pow(2) + (imag_mat.array().cast()).array().pow(2)).log10() * 10.0; + _flag = true; + } + else if (gdal_datatype == GDT_CInt16) { + Eigen::MatrixX> 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().array().pow(2) + complex_short_mat.imag().array().cast().pow(2)).log10() * 10.0; + _flag = true; + } + else if (gdal_datatype == GDT_CInt32) { + Eigen::MatrixX> 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().array().pow(2) + complex_short_mat.imag().array().cast().pow(2)).log10() * 10.0; + _flag = true; + } + else if (gdal_datatype == GDT_CFloat32) { + Eigen::MatrixX> 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().array().pow(2) + complex_short_mat.imag().array().cast().pow(2)).log10() * 10.0; + _flag = true; + } + else if (gdal_datatype == GDT_CFloat64) { + Eigen::MatrixX> 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().array().pow(2) + complex_short_mat.imag().array().cast().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 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 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 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 sign) +{ + return sqrt(pow(sign.real(),2)+pow(sign.imag(),2)); +} + +double LAMPTOOLAPI Complex2phase(std::complex 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().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); +} + + diff --git a/src/LAMPTool/SARBaseTool.h b/src/LAMPTool/SARBaseTool.h new file mode 100644 index 0000000..930c49d --- /dev/null +++ b/src/LAMPTool/SARBaseTool.h @@ -0,0 +1,66 @@ +#pragma once +/** +* 基础SARBaseTool 操作 +* +*/ +#include "LAMPToolAPI.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#define PI 3.141592653589793238462643383279 +#define PI_180 180/3.141592653589793238462643383279 +#define T180_PI 3.141592653589793238462643383279/180 +#define LIGHTSPEED 299792458 +#define LIGHESPEEDGHZ 0.299792458 + +#define Radians2Degrees(Radians) Radians*PI_180 +#define Degrees2Radians(Degrees) Degrees*T180_PI + + +// 判断是否需要输出为DLL +#define DLLOUT +int LAMPTOOLAPI Complex2dB(QString in_tiff, QString out_dB_path); +int LAMPTOOLAPI Amplitude2dB(QString in_tiff, QString out_dB_path); +Eigen::MatrixXd LAMPTOOLAPI Complex2dB(Eigen::MatrixXcd in_matrix); +Eigen::MatrixXd LAMPTOOLAPI Complex2dB(Eigen::MatrixXd in_matrix_real, Eigen::MatrixXd in_matrix_imag ); + +double LAMPTOOLAPI Complex2Amplitude(std::complex sign); +double LAMPTOOLAPI Complex2phase(std::complex sign);// 返回弧度制相位 + + + +Eigen::MatrixXd LAMPTOOLAPI Complex2Amplitude(Eigen::MatrixXcd in_matrix); +Eigen::MatrixXd LAMPTOOLAPI Complex2phase(Eigen::MatrixXcd in_matrix); +// ---------------------------------------------------------------------------------------------------------- +// 后向散射系数系统仿真模型 + +#ifndef DLLOUT + +#else + +#ifdef __cplusplus +extern "C" { +#endif +#define DllExport __declspec( dllexport ) + int __declspec(dllexport) Complex2dB_DLL(QString out_path, QString in_sar_path); + int __declspec(dllexport) Amplitude2dB_DLL(QString in_tiff, QString out_dB_path); +#ifdef __cplusplus +} +#endif + +#endif + diff --git a/src/LAMPTool/SARCalibration.cpp b/src/LAMPTool/SARCalibration.cpp new file mode 100644 index 0000000..92c0ade --- /dev/null +++ b/src/LAMPTool/SARCalibration.cpp @@ -0,0 +1,384 @@ +#include "SARCalibration.h" +#include "ImageOperatorBase.h" +#include "SARBaseTool.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +using namespace std; +using namespace Eigen; + +/** +* 输入数据是ENVI格式数据 +*/ + + + + + + +Eigen::Matrix2cd LAMPTOOLAPI CalibrationMatrix(Eigen::MatrixXcd in_matrix, double calibrationValue) +{ + return in_matrix.array() * calibrationValue; +} + +int LAMPTOOLAPI CalibrationSiglePolarSAR(QString out_path, QString in_sar_path, double calibrationValue) +{ + return CalibrationComplex(out_path, in_sar_path, calibrationValue); +} + +int LAMPTOOLAPI CalibrationComplex(const QString& out_path, const QString& in_sar_path, double calibrationValue) +{ + GDALAllRegister(); + std::shared_ptr rasterDataset = OpenDataset(in_sar_path); + int width = rasterDataset->GetRasterXSize(); + int height = rasterDataset->GetRasterYSize(); + int band_num = rasterDataset->GetRasterCount(); + double* gt = new double[6]; + std::shared_ptr gt_ptr(gt); + QString projectDef = rasterDataset->GetProjectionRef(); + GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType(); + bool _nGt_flag = false; + if (projectDef == "") { + _nGt_flag = true; + } + else { + _nGt_flag = false; + } + CreateDataset(out_path, height, width, 1, gt_ptr.get(), projectDef, GDT_CFloat64, _nGt_flag); + // 计算大小 + if (gdal_datatype == 0) { + return 1; + } + else if (gdal_datatype < 8) { + if (band_num != 2) { return 1; } + } + else if (gdal_datatype < 12) { + if (band_num != 1) { return 1; } + + } + else {} + int block_num = block_num_pre_memory(width, height, gdal_datatype, 2e9);// + block_num = block_num > height ? height : block_num; // 行数 + int line_num = block_num; + for (int i = 0; i < height; i= block_num+i) { + if (height - i < block_num) { + line_num = height - i; + } + else {} + // 构建矩阵块，使用eigen 进行矩阵计算，加速计算 + bool _flag = false; + Eigen::Matrix data_mat(line_num * width, 2);// 必须强制行优先 + if (gdal_datatype == GDT_Byte) { + Eigen::MatrixX real_mat(line_num*width,1); + Eigen::MatrixX 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() * calibrationValue).array(); + data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_UInt16) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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() * calibrationValue).array(); + data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_Int16) { + Eigen::MatrixX real_mat(line_num*width, 1); + Eigen::MatrixX 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() * calibrationValue).array(); + data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_UInt32) { + Eigen::MatrixX real_mat(line_num*width, 1); + Eigen::MatrixX 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() * calibrationValue).array(); + data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_Int32) { + Eigen::MatrixX real_mat(line_num*width, 1); + Eigen::MatrixX 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() * calibrationValue).array(); + data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + _flag = true; + } + //else if (gdal_datatype == GDT_UInt64) { + // Eigen::MatrixX real_mat(line_num*width, 1); + // Eigen::MatrixX 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() * calibrationValue).array(); + // data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + // _flag = true; + //} + //else if (gdal_datatype == GDT_Int64) { + // Eigen::MatrixX real_mat(line_num*width, 1); + // Eigen::MatrixX 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() * calibrationValue).array(); + // data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + // _flag = true; + //} + else if (gdal_datatype == GDT_Float32) { + Eigen::MatrixX real_mat(line_num*width, 1); + Eigen::MatrixX 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() * calibrationValue).array(); + data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_Float64) { + Eigen::MatrixX real_mat(line_num*width, 1); + Eigen::MatrixX 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() * calibrationValue).array(); + data_mat.col(1) =( imag_mat.array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_CInt16) { + Eigen::MatrixX> 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() * calibrationValue).array(); + data_mat.col(1) = (complex_short_mat.imag().array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_CInt32) { + Eigen::MatrixX> 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() * calibrationValue).array(); + data_mat.col(1) = (complex_short_mat.imag().array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_CFloat32) { + Eigen::MatrixX> 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() * calibrationValue).array(); + data_mat.col(1) = (complex_short_mat.imag().array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_CFloat64) { + Eigen::MatrixX> 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() * calibrationValue).array(); + data_mat.col(1) = (complex_short_mat.imag().array().cast() * calibrationValue).array(); + _flag = true; + } + else {} + // 保存数据 + if (_flag) { + // 定义赋值矩阵 + saveDataset(out_path, i,0,1, width, line_num, data_mat.data()); + } + else { + } + } + return -1; +} + +int LAMPTOOLAPI CalibrationComplex2dB(const QString& out_path, const QString& in_sar_path, double calibrationValue) +{ + GDALAllRegister(); + std::shared_ptr rasterDataset = OpenDataset(in_sar_path); + int width = rasterDataset->GetRasterXSize(); + int height = rasterDataset->GetRasterYSize(); + int band_num = rasterDataset->GetRasterCount(); + double* gt = new double[6]; + std::shared_ptr gt_ptr(gt); + QString projectDef = rasterDataset->GetProjectionRef(); + GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType(); + bool _nGt_flag = false; + if (projectDef == "") { + _nGt_flag = true; + } + else { + _nGt_flag = false; + } + CreateDataset(out_path, height, width, 1, gt_ptr.get(), projectDef, GDT_CFloat64, _nGt_flag); + // 计算大小 + if (gdal_datatype == 0) { + return 1; + } + else if (gdal_datatype < 8) { + if (band_num != 2) { return 1; } + } + else if (gdal_datatype < 12) { + if (band_num != 1) { return 1; } + + } + else {} + int block_num = block_num_pre_memory(width, height, gdal_datatype, 2e9);// + block_num = block_num > height ? height : block_num; // 行数 + int line_num = block_num; + for (int i = 0; i < height; i = block_num + i) { + if (height - i < block_num) { + line_num = height - i; + } + else {} + // 构建矩阵块，使用eigen 进行矩阵计算，加速计算 + bool _flag = false; + Eigen::Matrix data_mat(line_num * width, 2);// 必须强制行优先 + if (gdal_datatype == GDT_Byte) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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() * calibrationValue).array(); + data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_UInt16) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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() * calibrationValue).array(); + data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_Int16) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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() * calibrationValue).array(); + data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_UInt32) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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() * calibrationValue).array(); + data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_Int32) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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() * calibrationValue).array(); + data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + _flag = true; + } + //else if (gdal_datatype == GDT_UInt64) { + // Eigen::MatrixX real_mat(line_num * width, 1); + // Eigen::MatrixX 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() * calibrationValue).array(); + // data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + // _flag = true; + //} + //else if (gdal_datatype == GDT_Int64) { + // Eigen::MatrixX real_mat(line_num * width, 1); + // Eigen::MatrixX 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() * calibrationValue).array(); + // data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + // _flag = true; + //} + else if (gdal_datatype == GDT_Float32) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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() * calibrationValue).array(); + data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_Float64) { + Eigen::MatrixX real_mat(line_num * width, 1); + Eigen::MatrixX 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() * calibrationValue).array(); + data_mat.col(1) = (imag_mat.array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_CInt16) { + Eigen::MatrixX> 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() * calibrationValue).array(); + data_mat.col(1) = (complex_short_mat.imag().array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_CInt32) { + Eigen::MatrixX> 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() * calibrationValue).array(); + data_mat.col(1) = (complex_short_mat.imag().array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_CFloat32) { + Eigen::MatrixX> 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() * calibrationValue).array(); + data_mat.col(1) = (complex_short_mat.imag().array().cast() * calibrationValue).array(); + _flag = true; + } + else if (gdal_datatype == GDT_CFloat64) { + Eigen::MatrixX> 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() * calibrationValue).array(); + data_mat.col(1) = (complex_short_mat.imag().array().cast() * calibrationValue).array(); + _flag = true; + } + else {} + // 保存数据 + if (_flag) { + // 定义赋值矩阵 + + Eigen::Matrix data_out(line_num* width, 1);// 必须强制行优先 + data_out.col(0) = Complex2dB(data_mat.col(0).array(), data_mat.col(1).array()).array(); + saveDataset(out_path, i, 0, 1, width, line_num, data_out.data()); + } + else { + } + } + return -1; + + + + + return 0; +} + +int LAMPTOOLAPI CalibrationComplex_DLL(const QString& out_path, const QString& in_sar_path, double calibrationValue) +{ + return 0; +} + +int LAMPTOOLAPI CalibrationComplex2dB_DLL(const QString& out_path, const QString& in_sar_path, double calibrationValue) +{ + return 0; +} diff --git a/src/LAMPTool/SARCalibration.h b/src/LAMPTool/SARCalibration.h new file mode 100644 index 0000000..d51e273 --- /dev/null +++ b/src/LAMPTool/SARCalibration.h @@ -0,0 +1,52 @@ +#pragma once + + +#ifndef SARCALIBRATION_H +#define SARCALIBRATION_H +/** +* SAR定标函数 +* +*/ +#include "referenceHeader.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +using namespace std; +using namespace Eigen; +// 判断是否需要输出为DLL +#define DLLOUT + + +// 定标计算 +// 输入定标系统，保证数据能够正常定标 +Eigen::Matrix2cd LAMPTOOLAPI CalibrationMatrix(Eigen::MatrixXcd in_matrix, double calibrationValue); +int LAMPTOOLAPI CalibrationComplex(const QString& out_path, const QString& in_sar_path, double calibrationValue); +int LAMPTOOLAPI CalibrationComplex2dB(const QString& out_path, const QString& in_sar_path, double calibrationValue); +#ifndef DLLOUT + +#else + +#ifdef __cplusplus +extern "C" { +#endif + #define DllExport __declspec( dllexport ) + int __declspec(dllexport) CalibrationComplex_DLL(const QString& out_path, const QString& in_sar_path, double calibrationValue); + int __declspec(dllexport) CalibrationComplex2dB_DLL(const QString& out_path, const QString& in_sar_path, double calibrationValue); +#ifdef __cplusplus +} +#endif + + + +#endif + +#endif // !SARCALIBRATION_H \ No newline at end of file diff --git a/src/LAMPTool/SARCalibrationTool.cpp b/src/LAMPTool/SARCalibrationTool.cpp new file mode 100644 index 0000000..508e974 --- /dev/null +++ b/src/LAMPTool/SARCalibrationTool.cpp @@ -0,0 +1,87 @@ +// SARCalibrationTool.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。 +// +#include "LAMPToolAPI.h" +#include "SARCalibration.h" +#include "SARBaseTool.h" +#include +#include +#include + + + +int LAMPTOOLAPI SiglePolarCalibration(int argc, char* argv[]) { + qDebug() << "mode : 1 "<< "\n"; + QString in_tiff_path = ""; + QString out_path = ""; + double calibration = 0; + // params prase + for (int i = 0; i < argc; i++) { + if (string(argv[i]) == "-in") + { + in_tiff_path = argv[i + 1]; + qDebug() << "in file : " << in_tiff_path << "\n"; + } + else if (string(argv[i]) == "-out") { + out_path = argv[i + 1]; + qDebug() << "out file : " << in_tiff_path << "\n"; + } + else if (string(argv[i]) == "-calibrationConst") { + calibration = strtod(argv[i + 1], NULL); + qDebug() << "calibrationConst : " << calibration << "\n"; + } + else {} + } + // excute tool + return CalibrationComplex(out_path, in_tiff_path, calibration); +} + + +int LAMPTOOLAPI RasterComplex2dB(int argc, char* argv[]) { + qDebug() << "mode : 1 " << "\n"; + QString in_tiff_path = ""; + QString out_path = ""; + // params prase + for (int i = 0; i < argc; i++) { + if (string(argv[i]) == "-in") + { + in_tiff_path = argv[i + 1]; + qDebug() << "in file : " << in_tiff_path << "\n"; + } + else if (string(argv[i]) == "-out") { + out_path = argv[i + 1]; + qDebug() << "out file : " << in_tiff_path << "\n"; + } + else {} + } + // excute tool + return Complex2dB(in_tiff_path, out_path); +} + + +int LAMPTOOLAPI Test_SARCalibrationTool(int argc, char* argv[]) +{ + qDebug() << "calibration tool " << "\n"; + qDebug() << "sigle polsar Amg: SARCalibrationTool.exe 1 -in filepath -out filepath -calibrationConst 43" << "\n"; + qDebug() << "complex 2 dB : SARCalibrationTool.exe 1 -in filepath -out filepath" << "\n"; + if (argc == 1) { + qDebug() << "the number of params should be than 2" << "\n"; + return 2; + } + try { + int mode = atoi(argv[1]); + if (mode == 1) { + return SiglePolarCalibration(argc, argv); + } + else if (mode == 2) { + return RasterComplex2dB(argc, argv); + } + else {} + } + catch (std::exception ex) { + std::wcerr << ex.what() << "\n"; + } + + qDebug() << "Hello World!\n"; + return 0; +} + diff --git a/src/LAMPTool/SARImageBase.cpp b/src/LAMPTool/SARImageBase.cpp new file mode 100644 index 0000000..865c351 --- /dev/null +++ b/src/LAMPTool/SARImageBase.cpp @@ -0,0 +1,335 @@ +#include "LAMPToolAPI.h" +#include "ImageOperatorBase.h" +#include "SARBaseTool.h" +#include "SARImageBase.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "referenceHeader.h" +#include "OCCTBase.h" +#include "GeoOperator.h" +#include "SARImageBase.h" + + + + + +double LAMPTOOLAPI getRangeResolution(double startfreq, double endfreq) +{ + return LIGHTSPEED / 2 / std::abs(endfreq - startfreq); +} + +double LAMPTOOLAPI getAzimuthResolution(double AzAngleRange, double startFreq) +{ + return LIGHTSPEED / 2 / (AzAngleRange * startFreq); +} + +Eigen::MatrixXd LAMPTOOLAPI hammingWindows(size_t num) +{ + double alpha = 25. / 46; + double beta = 1 - alpha; + double scale = 1; + size_t m = num; + Eigen::MatrixXd x = Eigen::MatrixXd::Zero(m, 1); + for (int i = 0; i < m; i++) { + x(i, 0) = alpha - beta * std::cos(2 * PI * i / (num - 1)); + } + x = x.array() * scale; + return x; +} + + +Eigen::MatrixXd LAMPTOOLAPI Hanning(size_t Nf, size_t Nxa, double alpha) +{ + // n1 = double ? DOUBLE(n1In[0]) : FLOAT(n1In[0]) + double a = 2 * PI / Nf; // a = 2 * pi / N1 ;scale factor + double n2 = Nxa; // n2 = double ? DOUBLE(n2In[0]) : FLOAT(n2In[0]) + double b = 2 * pi / n2; // dim 2 scale fact + double one = 1.0; // double ? 1d : 1.0 + double temp_row = 0; + double temp_col = 0; + Eigen::MatrixXd result_arr = Eigen::MatrixXd::Zero(Nxa, Nf); // 结果 + for (int i = 0; i < Nxa; i++) { + for (int j = 0; j < Nf; j++) { + // index = double ? DINDGEN(n1) : FINDGEN(n1) ; Nf + temp_row = (alpha - one) * cos(j * a) + alpha; // (alpha-one) * cos(index*a) + alpha ;One row + // index = double ? DINDGEN(n2) : FINDGEN(n2) ; Nxa + temp_col = (alpha - one) * cos(i * b) + alpha; // col = (alpha-one) * cos(index*b) + alpha ;One column + result_arr(i, j) = temp_row * temp_col;// RETURN,(row # col) ;DINDGEN(n1)#DINDGEN(n2) + } + } + return result_arr; +} + + +Eigen::MatrixXcd LAMPTOOLAPI HammingWindows(Eigen::MatrixXcd FreqEcho) +{ + size_t data_cols = FreqEcho.cols(); + size_t data_rows = FreqEcho.rows(); + + // 进行进行hanmming windows + Eigen::MatrixXd hamming_Az = hammingWindows(data_cols); + Eigen::MatrixXd hamming_Rz = hammingWindows(data_rows); + for (int i = 0; i < data_cols; i++) { + for (int j = 0; j < data_rows; j++) { + FreqEcho(j, i) = FreqEcho(j, i) * hamming_Az(i, 0) * hamming_Rz(j, 0); + } + } + return FreqEcho; +} + +Eigen::MatrixXcd LAMPTOOLAPI InterpFreqEcho(Eigen::MatrixXcd freqEcho, Eigen::VectorXd SourceFreqlist, Eigen::VectorXd interpFreqlist) +{ + assert(freqEcho.cols() == SourceFreqlist.rows()); + + + + size_t data_rows = freqEcho.rows(); + size_t data_cols = freqEcho.cols(); + + size_t out_cols = interpFreqlist.rows(); + Eigen::MatrixXcd resultECHO = Eigen::MatrixXcd::Zero(data_rows, out_cols); + + +#pragma omp parallel for // NEW ADD + for (int i = 0; i < data_rows; i++) { // 单脉冲 + double* xs = new double[data_rows]; + for (int i = 0; i < data_cols; i++) { + xs[i] = SourceFreqlist(i); // 原始回波 + } + double* real_ys = new double[data_cols]; + double* imag_ys = new double[data_cols]; + + gsl_interp_accel* acc = gsl_interp_accel_alloc(); + const gsl_interp_type* t = gsl_interp_cspline_periodic; + gsl_spline* spline_real = gsl_spline_alloc(t, data_rows); + gsl_spline* spline_imag = gsl_spline_alloc(t, data_rows); + + for (int jj = 0; jj < data_cols; jj++) { + real_ys[jj] = freqEcho(jj, i).real(); + imag_ys[jj] = freqEcho(jj, i).imag(); + } + gsl_spline_init(spline_real, xs, real_ys, data_rows); + gsl_spline_init(spline_imag, xs, imag_ys, data_rows); + for (int j = 0; j < out_cols; j++) { // 受限于当前函数只能逐点插值，或者后期可修改为 gsl 直接插值 + double cx = interpFreqlist(j,0); + std::complex result( + gsl_spline_eval(spline_real, cx, acc),// 插值实部 + gsl_spline_eval(spline_imag, cx, acc)// 插值虚部 + ); + resultECHO(i, j) = result; + } + + gsl_spline_free(spline_real); + gsl_spline_free(spline_imag); + gsl_interp_accel_free(acc); + + delete[] xs; + delete[] real_ys, imag_ys; + } + return resultECHO; +} + +int LAMPTOOLAPI WriteComplexData2Amp_Arg(QString out_tiff_path, Eigen::MatrixXcd data) +{ + Eigen::MatrixXd amp_data = Eigen::MatrixXd::Zero(data.rows(), data.cols()); + Eigen::MatrixXd arg_data = Eigen::MatrixXd::Zero(data.rows(), data.cols()); + size_t row_count = data.rows(); + size_t col_count = data.cols(); + for (int i = 0; i < row_count; i++) { + for (int j = 0; j < col_count; j++) { + amp_data(i, j) = std::abs(data(i, j)); + arg_data(i, j) = std::arg(data(i, j)); + } + } + + + Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3); + + gdalImage image_tiff = CreategdalImage(out_tiff_path, row_count, col_count, 2, gt, "", false, true);// 注意这里保留仿真结果 + + image_tiff.saveImage(amp_data, 0, 0, 1); + image_tiff.saveImage(arg_data, 0, 0, 2); + + return 0; +} + +int LAMPTOOLAPI WriteComplexData2AmpdB_Arg(QString out_tiff_path, Eigen::MatrixXcd data) +{ + Eigen::MatrixXd amp_data = Eigen::MatrixXd::Zero(data.rows(), data.cols()); + Eigen::MatrixXd arg_data = Eigen::MatrixXd::Zero(data.rows(), data.cols()); + size_t row_count = data.rows(); + size_t col_count = data.cols(); + for (int i = 0; i < row_count; i++) { + for (int j = 0; j < col_count; j++) { + amp_data(i, j) = 10.0*std::log10(std::abs(data(i, j))); + arg_data(i, j) = std::arg(data(i, j)); + } + } + + + Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3); + + gdalImage image_tiff = CreategdalImage(out_tiff_path, row_count, col_count, 2, gt, "", false, true);// 注意这里保留仿真结果 + + image_tiff.saveImage(amp_data, 0, 0, 1); + image_tiff.saveImage(arg_data, 0, 0, 2); + + return 0; +} + +size_t LAMPTOOLAPI nextpow2(size_t num) +{ + double n = std::round(std::log10(num * 1.0) / std::log10(2.0)); + return pow(2, (size_t)std::ceil(n) + 1); +} + + +Eigen::MatrixXcd LAMPTOOLAPI IFFTW1D(Eigen::MatrixXcd ECHOdata) +{ + size_t data_rows = ECHOdata.rows(); + size_t data_cols = ECHOdata.cols(); + // 频域 转换到 时域 + size_t n2 = nextpow2(data_cols) ; + //qDebug() << data_rows << "," << data_cols << "," << n2 << "\n"; + // 补零 + Eigen::MatrixXcd ExpandECHO = Eigen::MatrixXcd::Zero(data_rows, n2); + for (int i = 0; i < data_rows; i++) { + for (int j = 0; j < data_cols; j++) { + ExpandECHO(i, j) = ECHOdata(i, j); + } + } + // 结果 + Eigen::MatrixXcd echoTime = Eigen::MatrixXcd::Zero(data_rows, n2); + + fftw_complex* din = (fftw_complex*)fftw_malloc(sizeof(double) * n2 * 2); + fftw_complex* dout = (fftw_complex*)fftw_malloc(sizeof(double) * n2 * 2); + fftw_plan p; + double n = 1.0 / n2; // 因为 fftw 的逆傅里叶变换，没有除于N,所以实际上是原来的N倍 + //逐行进行fftw + p = fftw_plan_dft_1d(n2, din, dout, FFTW_BACKWARD, FFTW_MEASURE);//FFTW_BACKWARD 逆， + fftw_execute(p); //FFTW_MEASURE 估计最优变换方法，后面复用变换方案 + // 傅里叶 计算进度 + for (int i = 0; i < data_rows; i++) { + for (int j = 0; j < n2; j++) { // 初始化输入 + din[j][0] = std::real(ExpandECHO(i, j)); + din[j][1] = std::imag(ExpandECHO(i, j)); + } + // 构建fftw 任务 + fftw_execute(p); + // 保存结果 + for (int j = 0; j < n2; j++) { //读取结果 + echoTime(i, j) = std::complex(dout[j][0], dout[j][1]) * n; // 每次变换后都除于n + } + printf("\rIFFT[%.2lf%%]...:", i * 100.0 / (data_rows - 1)); + } + fftw_destroy_plan(p); + fftw_free(din); + fftw_free(dout); + qDebug() << "\n"; + qDebug() << "IFFT over" << "\n"; + return echoTime; +} + +Eigen::MatrixXcd LAMPTOOLAPI FFTW1D(Eigen::MatrixXcd ECHO) +{ + size_t data_rows = ECHO.rows(); + size_t data_cols = ECHO.cols(); + // 频域 转换到 时域 + size_t n2 = data_cols; + // 补零 + Eigen::MatrixXcd ExpandECHO = Eigen::MatrixXcd::Zero(data_rows, n2); + for (int i = 0; i < data_rows; i++) { + for (int j = 0; j < data_cols; j++) { + ExpandECHO(i, j) = ECHO(i, j); + } + } + // 结果 + Eigen::MatrixXcd echofreq = Eigen::MatrixXcd::Zero(data_rows, n2); + fftw_complex* din = (fftw_complex*)fftw_malloc(sizeof(double) * n2 * 2); + fftw_complex* dout = (fftw_complex*)fftw_malloc(sizeof(double) * n2 * 2); + fftw_plan p; + //逐行进行fftw + p = fftw_plan_dft_1d(n2, din, dout, FFTW_FORWARD, FFTW_MEASURE);//FFTW_FORWARD 顺， + fftw_execute(p); //FFTW_MEASURE 估计最优变换方法，后面复用变换方案 + for (int i = 0; i < data_rows; i++) { + for (int j = 0; j < n2; j++) { // 初始化输入 + din[j][0] = std::real(ExpandECHO(i, j)); + din[j][1] = std::imag(ExpandECHO(i, j)); + } + // 构建fftw 任务 + fftw_execute(p); + // 保存结果 + for (int j = 0; j < n2; j++) { //读取结果 + echofreq(i, j) = std::complex(dout[j][0], dout[j][1]); + } + printf("\rFFT[%.2lf%%]...:", i * 100.0 / (data_rows - 1)); + } + fftw_destroy_plan(p); + fftw_free(din); + fftw_free(dout); + qDebug() << "\n"; + qDebug() << "FFT over" << "\n"; + return echofreq; +} + + + +Eigen::MatrixXcd LAMPTOOLAPI FFTW2D(Eigen::MatrixXcd ECHO) +{ + // 获取数据的尺寸 + int rows = ECHO.rows(); + int cols = ECHO.cols(); + + // 创建FFTW输入和输出数组 + fftw_complex* in = reinterpret_cast(ECHO.data()); + fftw_complex* out = static_cast(fftw_malloc(sizeof(fftw_complex) * rows * cols)); + + // 创建FFTW计划，执行傅立叶变换 + fftw_plan plan = fftw_plan_dft_2d(rows, cols, in, out, FFTW_FORWARD, FFTW_ESTIMATE); + fftw_execute(plan); + + // 打印傅立叶变换结果 + Eigen::MatrixXcd result(rows, cols); + result = Map(reinterpret_cast*>(out), rows, cols); + //qDebug() << "Fourier Transform Result: " << "\n" << result << "\n"; + + // 销毁FFTW计划和内存 + fftw_destroy_plan(plan); + fftw_free(out); + qDebug() << "\n"; + qDebug() << "FFT over" << "\n"; + return result; +} + + +Eigen::MatrixXcd LAMPTOOLAPI fftshift(Eigen::MatrixXcd X) { + int rows = X.rows(); + int cols = X.cols(); + int rows_half = rows / 2; + int cols_half = cols / 2; + + Eigen::MatrixXcd tmp = X.topLeftCorner(rows_half, cols_half); + X.topLeftCorner(rows_half, cols_half) = X.bottomRightCorner(rows_half, cols_half); + X.bottomRightCorner(rows_half, cols_half) = tmp; + + tmp = X.topRightCorner(rows_half, cols_half); + X.topRightCorner(rows_half, cols_half) = X.bottomLeftCorner(rows_half, cols_half); + X.bottomLeftCorner(rows_half, cols_half) = tmp; + + return X; +} \ No newline at end of file diff --git a/src/LAMPTool/SARImageBase.h b/src/LAMPTool/SARImageBase.h new file mode 100644 index 0000000..2cbb298 --- /dev/null +++ b/src/LAMPTool/SARImageBase.h @@ -0,0 +1,95 @@ +#pragma once +#ifndef SARIMAGEBASE_H +#define SARIMAGEBASE_H +/** +* 成像过程中常用的计算函数 +****/ +#include "referenceHeader.h" +#include "SARBaseTool.h" +// ------------------------------------------- 成像 公用库---------------------------------------------------------- +// SAR 常用计算的方法 +double LAMPTOOLAPI getRangeResolution(double startfreq, double endfreq); + +/// +/// 方位向分辨率计算 +/// +/// 弧度制，方位角变化范围（>=0） +/// 起始频率 +/// +double LAMPTOOLAPI getAzimuthResolution(double AzAngleRange, double startFreq); +// 根据输入数量构建汉明窗 +Eigen::MatrixXd LAMPTOOLAPI hammingWindows(size_t num); + + +// 频率域汉明窗口 + +/// +/// 频率域汉明窗口 +/// +/// 频率点数 +/// 脉冲数量 +/// 权重 +/// 汉明窗权重 +Eigen::MatrixXd LAMPTOOLAPI Hanning(size_t Nf, size_t Nxa, double alpha = 0.54); + +/// +/// 获取2的平方数,大于2 +/// +/// +/// +size_t LAMPTOOLAPI nextpow2(size_t num); + +/// +/// 对回波沿着行进行一维傅里叶逆变换 +/// +/// +/// +Eigen::MatrixXcd LAMPTOOLAPI IFFTW1D(Eigen::MatrixXcd ECHO); + + +/// +/// 对回波沿着行进行一维傅里叶变换 +/// +/// +/// +Eigen::MatrixXcd LAMPTOOLAPI FFTW1D(Eigen::MatrixXcd echo); +Eigen::MatrixXcd LAMPTOOLAPI FFTW2D(Eigen::MatrixXcd ECHO); + +Eigen::MatrixXcd LAMPTOOLAPI fftshift(Eigen::MatrixXcd X); +/// +/// 对回波 加 hamming 窗 +/// +/// 频域的回波数据（PRFnum,Freqnum) +/// +Eigen::MatrixXcd LAMPTOOLAPI HammingWindows(Eigen::MatrixXcd FreqEcho); + + +/// +/// 对回波进行一维插值(并行），不能用来处理极坐标 +/// 插值方式为 实部虚部分别进行线性插值 +/// +/// 回波数据,（行数，列数）:（PRFNUM,freqNUM） +/// 回波的 原始频率列表 +/// 回波的 插值频率列表 +/// +Eigen::MatrixXcd LAMPTOOLAPI InterpFreqEcho(Eigen::MatrixXcd freqEcho, Eigen::VectorXd SourceFreqlist, Eigen::VectorXd interpFreqlist); + + +/// +/// 保存复数矩阵数据 为tiff， 其中，band 1: amp (linear) angle(radia) +/// +/// +/// +/// +int LAMPTOOLAPI WriteComplexData2Amp_Arg(QString out_tiff_path, Eigen::MatrixXcd data); + + +/// +/// 保存复数矩阵数据 为tiff， 其中，band 1: amp (dB) angle(radia) +/// +/// +/// +/// +int LAMPTOOLAPI WriteComplexData2AmpdB_Arg(QString out_tiff_path, Eigen::MatrixXcd data); +// ---------------------------------------------------------------------------------------------------------- +#endif diff --git a/src/LAMPTool/interpolation.cpp b/src/LAMPTool/interpolation.cpp new file mode 100644 index 0000000..bbf1b84 --- /dev/null +++ b/src/LAMPTool/interpolation.cpp @@ -0,0 +1,13 @@ +#include "interpolation.h" +#include +#include +#include +#include +#include +#include + + + +namespace LampInterpolation { + +} \ No newline at end of file diff --git a/src/LAMPTool/interpolation.h b/src/LAMPTool/interpolation.h new file mode 100644 index 0000000..a93a8ab --- /dev/null +++ b/src/LAMPTool/interpolation.h @@ -0,0 +1,69 @@ +#pragma once + +#ifndef INTERPOLATION_H +#define INTERPOLATION_H + +#include +#include +#include +#include +#include +#include +#include + + +namespace LampInterpolation { + enum interpolationtype { + nearest, + linear, + cubic, + }; + + + /// + /// 矩阵为 1xn + /// + template + std::complex interpolation(Eigen::MatrixX>& echo, double& index, interpolationtype methodtype) + { + assert(echo.rows() != 1 || echo.cols() >= index || index < 0); // 断言 + if (methodtype == interpolationtype::linear) { + return interpolationLinear(echo, index); + } + else if (methodtype == interpolationtype::cubic) { + + } + else if (methodtype == interpolationtype::nearest) { + return interpolationNearest(echo, index); + } + else {} + return std::complex(0, 0); + }; + + template + std::complex interpolationLinear(Eigen::MatrixX>& echo, double& index) + { + size_t last_ids = size_t(std::floor(index)); + size_t next_ids = size_t(std::ceil(index)); + + std::complex last_value = echo(1, last_ids); + std::complex next_value = echo(1, next_ids); + + // 实部，虚部同时插值 + double real = last_value.real() + ((next_value.real() - last_value.real()) / (next_ids - last_ids)) * (index - last_ids); + double imag = last_value.imag() + ((next_value.imag() - last_value.imag()) / (next_ids - last_ids)) * (index - last_ids); + + return std::complex(T(real), T(imag)); + }; + + template + std::complex interpolationNearest(Eigen::MatrixX>& echo, double& index) + { + size_t nearest_ids = size_t(std::round(index)); + return echo(1, nearest_ids); + }; + +} + + +#endif \ No newline at end of file diff --git a/src/PluginWBFZExchangePlugin/CMakeLists.txt b/src/PluginWBFZExchangePlugin/CMakeLists.txt new file mode 100644 index 0000000..40fe298 --- /dev/null +++ b/src/PluginWBFZExchangePlugin/CMakeLists.txt @@ -0,0 +1,153 @@ +#----------------------------------------------------------------------------- +# 将点云与 mesh 表达合并，不再单独处理点云 +# 加载： 点云-> mesh +# 处理： mesh->点云-> 处理模块 -> mesh +# +#----------------------------------------------------------------------------- + +#----------------------------------------------------------------------------- +# 头文件搜索路径 +#----------------------------------------------------------------------------- +include_directories(${CMAKE_CURRENT_SOURCE_DIR}/..) + +# boost +include_directories(D:/vcpkg/installed/x64-windows/include) + +# pcl +include_directories(C:/PCL/3rdParty/FLANN/include) +include_directories(C:/PCL/3rdParty/VTK/include/vtk-9.3) +include_directories(C:/PCL/include/pcl-1.14) + +# qt5 +include_directories(C:/Qt/5.15.2/msvc2019_64/include/QtQml) + +#----------------------------------------------------------------------------- +# 链接库 +#----------------------------------------------------------------------------- + +# pcl +link_directories("C:/PCL/3rdParty/FLANN/lib") +link_directories("C:/VTK/lib") +link_directories("C:/PCL/lib") + + +#----------------------------------------------------------------------------- +# 自动添加include目录 +#----------------------------------------------------------------------------- +set(CMAKE_INCLUDE_CURRENT_DIR ON) + + + + +#----------------------------------------------------------------------------- +# 添加资源文件 +#----------------------------------------------------------------------------- +set(_qrc "${CMAKE_CURRENT_SOURCE_DIR}/../qrc/qianfan.qrc") +set(_lang "${CMAKE_CURRENT_SOURCE_DIR}/../qrc/translations.qrc") +qt5_add_resources(_resource ${_qrc} ${_lang}) + +#----------------------------------------------------------------------------- +# 源码扫描 +#----------------------------------------------------------------------------- +file(GLOB _ui "*.ui") +file(GLOB _header "*.h") +file(GLOB _source "*.cpp") +qt5_wrap_ui(_interface ${_ui}) + +#----------------------------------------------------------------------------- +# 添加动态库目标 +#----------------------------------------------------------------------------- +add_library(PluginWBFZExchangePlugin + ${_resource} + ${_interface} + ${_header} + ${_source} +) + + +#----------------------------------------------------------------------------- +# 链接库 +#----------------------------------------------------------------------------- + +# pcl +link_directories("C:/PCL/3rdParty/FLANN/lib") +link_directories("C:/VTK/lib") +link_directories("C:/PCL/lib") + +find_package(Qt5 REQUIRED COMPONENTS Core Quick Sql Core Xml Opengl Gui Svg Xmlpatterns Uitools Widgets Qml Printsupport Sensors Quickwidgets Quick Concurrent Openglextensions Charts Datavisualization) +find_package(PCL ) + +# boost +include_directories(D:/vcpkg/installed/x64-windows/include) + +# pcl +include_directories(C:/PCL/3rdParty/FLANN/include) +include_directories(C:/PCL/3rdParty/VTK/include/vtk-9.3) +include_directories(C:/PCL/include/pcl-1.14) + + +include_directories(${PCL_INCLUDE_DIRS}) +link_directories(${PCL_LIBRARY_DIRS}) +add_definitions(${PCL_DEFINITIONS}) + + + +#----------------------------------------------------------------------------- +# 添加接口声明宏 +#----------------------------------------------------------------------------- +target_compile_definitions(PluginWBFZExchangePlugin PRIVATE "NOMINMAX") # 禁用vc++ 中的min max 定义，避免与 pcl 的中 std::min，pcl::max 函数冲突 +target_compile_definitions(PluginWBFZExchangePlugin PRIVATE "WBFZ_API") + +list(APPEND _depend_library + LAMPTool WBCLFZSystemModule Common PointCloudOperator PythonModule DataProperty MeshData Material Geometry BCBase ConfigOptions SelfDefObject ModelData ModuleBase PluginManager GeometryCommand GeometryWidgets IO MainWidgets MainWindow +) + + + + +list(APPEND _runtimes_libraries + LAMPCAE::CGNS Qt5::Core Qt5::Gui Qt5::Widgets Qt5::Xml +) + +# 添加 点云 pcl +list(APPEND _runtimes_libraries + +) + + +# 添加 VTK +list(APPEND _runtimes_libraries + VTK::CommonColor VTK::CommonComputationalGeometry VTK::CommonCore VTK::CommonDataModel VTK::CommonExecutionModel VTK::CommonMath VTK::CommonMisc VTK::CommonSystem VTK::CommonTransforms VTK::DICOMParser VTK::FiltersCore VTK::FiltersGeneral VTK::FiltersGeometry VTK::FiltersHybrid VTK::FiltersModeling VTK::FiltersSources VTK::IOCore VTK::IOChemistry VTK::IOGeometry VTK::IOImage VTK::IOLegacy VTK::ImagingCore VTK::ImagingSources VTK::RenderingCore VTK::RenderingLOD VTK::doubleconversion VTK::jpeg VTK::jsoncpp VTK::lz4 VTK::lzma VTK::metaio VTK::png VTK::pugixml VTK::sys VTK::tiff VTK::zlib +) + +#----------------------------------------------------------------------------- +# 链接依赖库 +#----------------------------------------------------------------------------- +target_link_libraries(PluginWBFZExchangePlugin PRIVATE + ${_runtimes_libraries} + ${_depend_library} + ${PCL_LIBRARIES} +) + +#----------------------------------------------------------------------------- +# 添加依赖关系 +#----------------------------------------------------------------------------- +add_dependencies(PluginWBFZExchangePlugin ${_depend_library}) + +#----------------------------------------------------------------------------- +# 添加运行时依赖关系 +#----------------------------------------------------------------------------- +set(LAMPCAE_PluginWBFZExchangePlugin_Runtimes_Libraries ${_runtimes_libraries} PARENT_SCOPE) + +#----------------------------------------------------------------------------- +# 设置插件的输出目录 +#----------------------------------------------------------------------------- +set_target_properties(PluginWBFZExchangePlugin + PROPERTIES + ARCHIVE_OUTPUT_DIRECTORY_RELEASE $/plugins + LIBRARY_OUTPUT_DIRECTORY_RELEASE $/plugins + RUNTIME_OUTPUT_DIRECTORY_RELEASE $/plugins + ARCHIVE_OUTPUT_DIRECTORY_DEBUG $/plugins + LIBRARY_OUTPUT_DIRECTORY_DEBUG $/plugins + RUNTIME_OUTPUT_DIRECTORY_DEBUG $/plugins +) \ No newline at end of file diff --git a/src/PluginWBFZExchangePlugin/PointCloudThreadBase.cpp b/src/PluginWBFZExchangePlugin/PointCloudThreadBase.cpp new file mode 100644 index 0000000..fd5cb29 --- /dev/null +++ b/src/PluginWBFZExchangePlugin/PointCloudThreadBase.cpp @@ -0,0 +1,86 @@ +#include "PointCloudThreadBase.h" +#include "MainWindow/MainWindow.h" +#include "MainWindow/SubWindowManager.h" +#include "PythonModule/PyAgent.h" +#include "MeshData/meshSingleton.h" +#include "MeshData/meshKernal.h" +#include "Common/Types.h" +#include "Common/DebugLogger.h" + +namespace WBFZ +{ + PointCloudThreadBase::PointCloudThreadBase(const QString &fileName, PointCloudOperation operation, GUI::MainWindow *mw) : ModuleBase::ThreadTask(mw), + _fileName(fileName), + _operation(operation) + { + } + + PointCloudThreadBase::~PointCloudThreadBase() + { + } + + void PointCloudThreadBase::defaultMeshFinished() + { + + if (_threadRuning) + { + QString information{}; + ModuleBase::Message msg; + if (_operation == POINTCLOUD_FILTER) + { + if (_success) + { + emit _mainwindow->updateMeshTreeSig(); + emit _mainwindow->updateSetTreeSig(); + emit _mainwindow->updateActionStatesSig(); + // emit _mainwindow->updateActionsStatesSig(); + emit _mainwindow->getSubWindowManager()->openPreWindowSig(); + emit _mainwindow->updatePreMeshActorSig(); + information = QString("Successful Import Mesh From \"%1\"").arg(_fileName); + msg.type = Common::Message::Normal; + msg.message = information; + auto meshdata = MeshData::MeshData::getInstance(); + // meshdata->generateDisplayDataSet(); + const int nk = meshdata->getKernalCount(); + if (nk <= 0) + return; + auto k = meshdata->getKernalAt(nk - 1); + if (k != nullptr) + k->setPath(_fileName); + } + else + { + information = QString("Failed Filter From \"%1\"").arg(_fileName); + msg.type = Common::Message::Error; + msg.message = information; + } + } + else if (_operation == POINTCLOUD_MESH) + { + if (_success) + { + information = QString("Successful resurface Mesh To \"%1\"").arg(_fileName); + msg.type = Common::Message::Normal; + msg.message = information; + } + else + { + information = QString("Failed resurface Mesh To \"%1\"").arg(_fileName); + msg.type = Common::Message::Error; + msg.message = information; + } + } + emit showInformation(information); + emit _mainwindow->printMessageToMessageWindow(msg); + } + // emit showInformation(information); + // emit _mainwindow->printMessageToMessageWindow(msg); + ModuleBase::ThreadTask::threadTaskFinished(); + Py::PythonAgent::getInstance()->unLock(); + } + void PointCloudThreadBase::setThreadRunState(bool flag){ + _success= flag; + + } + +} \ No newline at end of file diff --git a/src/PluginWBFZExchangePlugin/PointCloudThreadBase.h b/src/PluginWBFZExchangePlugin/PointCloudThreadBase.h new file mode 100644 index 0000000..512dddf --- /dev/null +++ b/src/PluginWBFZExchangePlugin/PointCloudThreadBase.h @@ -0,0 +1,33 @@ +#ifndef _MESHTHREADBASE_H_ +#define _MESHTHREADBASE_H_ + +#include "ModuleBase/ThreadTask.h" +#include "WBFZExchangePluginAPI.h" +#include "WBFZExchangePlugin.h" + +namespace GUI +{ + class MainWindow; +} + +namespace WBFZ +{ + class WBFZAPI PointCloudThreadBase : public ModuleBase::ThreadTask + { + public: + PointCloudThreadBase(const QString &fileName, PointCloudOperation operation, GUI::MainWindow *mw); + virtual ~PointCloudThreadBase(); + + virtual void run() = 0; + void defaultMeshFinished(); + void setThreadRunState(bool); + + + private: + bool _success{false}; + QString _fileName; + PointCloudOperation _operation; + }; +} + +#endif \ No newline at end of file diff --git a/src/PluginWBFZExchangePlugin/WBFZExchangePlugin.cpp b/src/PluginWBFZExchangePlugin/WBFZExchangePlugin.cpp new file mode 100644 index 0000000..9e94d8c --- /dev/null +++ b/src/PluginWBFZExchangePlugin/WBFZExchangePlugin.cpp @@ -0,0 +1,55 @@ + +#include "IO/IOConfig.h" + +#include "ModelData/modelDataBase.h" +#include "ModuleBase/ThreadControl.h" +#include "ModuleBase/ThreadTaskManager.h" +#include "WBFZExchangePlugin.h" + +#include + + + +namespace WBFZ { + GUI::MainWindow* WBFZ::WBFZExchangePlugin::_mainwindow = nullptr; + + WBFZExchangePlugin::WBFZExchangePlugin(GUI::MainWindow* m) + { + _describe = "WBFZExchangePlugin Installed Successfully"; + _mainwindow = m; + } + + bool WBFZExchangePlugin::install() + { + IO::IOConfigure::RegisterMeshImporter("CGNS(*.cgns)", CGNSimportMesh); + +// IO::IOConfigure::RegisterMeshExporter("CGNS(*.cgns)", CGNSexportMesh); + + + + return true; + } + + bool WBFZExchangePlugin::uninstall() + { + IO::IOConfigure::RemoveMeshImporter("CGNS(*.cgns)"); + + IO::IOConfigure::RemoveMeshExporter("CGNS(*.cgns)"); + + return true; + } + + void WBFZExchangePlugin::reTranslate(QString) {} + + GUI::MainWindow* WBFZExchangePlugin::getMWpt() + { + return _mainwindow; + } +} // namespace WBFZ + +void Register(GUI::MainWindow* m, QList* ps) +{ + Plugins::PluginBase* p = new WBFZ::WBFZExchangePlugin(m); + ps->append(p); +} + diff --git a/src/PluginWBFZExchangePlugin/WBFZExchangePlugin.h b/src/PluginWBFZExchangePlugin/WBFZExchangePlugin.h new file mode 100644 index 0000000..e00079a --- /dev/null +++ b/src/PluginWBFZExchangePlugin/WBFZExchangePlugin.h @@ -0,0 +1,42 @@ +#ifndef _WBFZEXCHANGEPLUGIN_H_ +#define _WBFZEXCHANGEPLUGIN_H_ + +#include "WBFZExchangePluginAPI.h" +#include "PluginManager/pluginBase.h" + +namespace WBFZ +{ + enum PointCloudOperation + { + POINTCLOUD_NONE, + POINTCLOUD_FILTER, + POINTCLOUD_MESH + }; + + class WBFZAPI WBFZExchangePlugin : public Plugins::PluginBase + { + public: + WBFZExchangePlugin(GUI::MainWindow* m); + ~WBFZExchangePlugin() = default; + bool install(); + bool uninstall(); + void reTranslate(QString); + static GUI::MainWindow* getMWpt(); + + private: + static GUI::MainWindow* _mainwindow; + }; +} + +extern "C" +{ + void WBFZAPI Register(GUI::MainWindow* m, QList* plugs); + //函数返回值是无效的,不要通过返回值判断 + bool WBFZAPI CGNSimportMesh(QString AbFileName, int modelId); + bool WBFZAPI MSHimportMesh(QString AbFileName, int modelId); + + + +} + +#endif \ No newline at end of file diff --git a/src/PluginWBFZExchangePlugin/WBFZExchangePluginAPI.h b/src/PluginWBFZExchangePlugin/WBFZExchangePluginAPI.h new file mode 100644 index 0000000..7035206 --- /dev/null +++ b/src/PluginWBFZExchangePlugin/WBFZExchangePluginAPI.h @@ -0,0 +1,13 @@ +#ifndef _WBFZAPI_H__ +#define _WBFZAPI_H__ + +#include + + +#if defined(WBFZ_API) +#define WBFZAPI Q_DECL_EXPORT +#else +#define WBFZAPI Q_DECL_IMPORT +#endif + +#endif