diff --git a/.gitignore b/.gitignore index 2a16f91..209932b 100644 --- a/.gitignore +++ b/.gitignore @@ -9,7 +9,11 @@ CMakeUserPresets.json /vcpkg/ # ---- external source-built deps (方案②-修订: VTK 源码编到 install 前缀) ---- -/external/ +# 用 /external/* 而非 /external/ 忽略目录内容，使下方 vendored 子目录可被例外重新纳入 +# (git 无法重新纳入被整体忽略目录内的文件)。 +/external/* +# 例外：vendored 3DGPRViewer 数据生成链(原样拷贝的算法源码，版权自有)需入库。 +!/external/gpr3dviewer/ # ---- Visual Studio / IDE ---- .vs/ diff --git a/CMakeLists.txt b/CMakeLists.txt index 57c0e2f..ac9ea1a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -76,10 +76,17 @@ if(EXISTS "${CMAKE_SOURCE_DIR}/external/qwt-src/src") include("${CMAKE_SOURCE_DIR}/cmake/qwt.cmake") endif() +# vendored 3DGPRViewer 数据生成链（原样拷贝，算法零改动）：geopro_gpr3dv 静态库。 +# 链：多通道 .iprh/.iprb → GPRDataModel 立方体 → RadarProcessor 处理。生产管线 A 地基。 +add_subdirectory(external/gpr3dviewer) + add_subdirectory(src) # POC-B headless 度量 CLI（gpr_poc）。链 io_gpr/core/store/render，在真实数据上跑端到端度量。 add_subdirectory(tools/gpr_poc) +# gpr3dv 冒烟 CLI：走 vendored 原版 API（loadImpulseMultiChannel → buildVolumeData → runPipeline）。 +add_subdirectory(tools/gpr3dv_smoke) + enable_testing() add_subdirectory(tests) diff --git a/external/gpr3dviewer/CMakeLists.txt b/external/gpr3dviewer/CMakeLists.txt new file mode 100644 index 0000000..140e85d --- /dev/null +++ b/external/gpr3dviewer/CMakeLists.txt @@ -0,0 +1,49 @@ +# ===================================================================== +# geopro_gpr3dv —— vendored 3DGPRViewer 数据生成链（原样拷贝，算法零改动） +# +# 链路：多通道 _Axx.iprh/.iprb +# → IprhParser::loadImpulseMultiChannel → GPRDataModel（traces） +# → GPRDataModel::buildVolumeData → volumeData[ch][trace][sample] +# → RadarProcessor::runPipeline → 处理后 GPRDataModel +# +# 版权属用户（lanbingtech 自有），可拷入本仓库 vendored 隔离。 +# 仅链 Qt6::Core（QVector/QString/QFile/QVector3D/QJson）+ OpenMP + kissfft。 +# 不含 UI/Sql/Network（这些在原版属 GPRWidget/mainwindow，未拷）。 +# ===================================================================== + +# 根工程仅 enable CXX；kissfft 是 C 源（kiss_fft.c/kiss_fftr.c），需显式启用 C 语言， +# 否则 .c 文件被 CMake 忽略，导致 kiss_fftr/kiss_fftr_alloc 链接缺失。 +enable_language(C) + +find_package(OpenMP) + +add_library(geopro_gpr3dv STATIC + IprhParser.cpp + ImpulseMultiChannelConverter.cpp + Rd3Parser.cpp + RadarProcessor.cpp + PerformanceLogger.cpp + third_party/kissfft/kiss_fft.c + third_party/kissfft/kiss_fftr.c +) + +# 头文件目录：本目录（GPRDataModel.h 等）+ kissfft（RadarProcessor.cpp 内 #include "kiss_fftr.h"）。 +target_include_directories(geopro_gpr3dv PUBLIC + ${CMAKE_CURRENT_SOURCE_DIR} + ${CMAKE_CURRENT_SOURCE_DIR}/third_party/kissfft +) + +# Qt6::Gui 为 QVector3D（GPRDataModel.h 的道三维位置）所需；Core 提供 QVector/QString/QFile/QJson。 +target_link_libraries(geopro_gpr3dv PUBLIC Qt6::Core Qt6::Gui) + +if(OpenMP_CXX_FOUND) + target_link_libraries(geopro_gpr3dv PUBLIC OpenMP::OpenMP_CXX) +endif() + +target_compile_features(geopro_gpr3dv PRIVATE cxx_std_17) + +# vendored 第三方代码：关 Qt 自动工具（无 QObject/moc 需求），并放宽告警（不改算法、不为告警动代码）。 +set_target_properties(geopro_gpr3dv PROPERTIES AUTOMOC OFF AUTOUIC OFF AUTORCC OFF) +if(MSVC) + target_compile_options(geopro_gpr3dv PRIVATE /W0) +endif() diff --git a/external/gpr3dviewer/GPRDataModel.h b/external/gpr3dviewer/GPRDataModel.h new file mode 100644 index 0000000..664eb16 --- /dev/null +++ b/external/gpr3dviewer/GPRDataModel.h @@ -0,0 +1,297 @@ +#ifndef GPRDATAMODEL_H +#define GPRDATAMODEL_H + +#include +#include +#include +#include +#include +#include +#include "SurveyGeometry.h" + +struct RadarTrace { + QVector amplitudes; // 【重要】改为 short (16-bit)，之前是 float + double startTime = 0.0; + double timeStep = 0.0; + // 三维位置信息 + QVector3D position; // 该道的三维空间位置 + int channelNumber = 0; // 通道号 + + RadarTrace() = default; + RadarTrace(const RadarTrace& other) + : amplitudes(QVector(other.amplitudes.cbegin(), other.amplitudes.cend())) + , startTime(other.startTime) + , timeStep(other.timeStep) + , position(other.position) + , channelNumber(other.channelNumber) + {} + RadarTrace& operator=(const RadarTrace& other) { + if (this != &other) { + amplitudes = QVector(other.amplitudes.cbegin(), other.amplitudes.cend()); + startTime = other.startTime; + timeStep = other.timeStep; + position = other.position; + channelNumber = other.channelNumber; + } + return *this; + } + RadarTrace(RadarTrace&& other) noexcept + : amplitudes(std::move(other.amplitudes)) + , startTime(other.startTime) + , timeStep(other.timeStep) + , position(other.position) + , channelNumber(other.channelNumber) + {} + RadarTrace& operator=(RadarTrace&& other) noexcept { + if (this != &other) { + amplitudes = std::move(other.amplitudes); + startTime = other.startTime; + timeStep = other.timeStep; + position = other.position; + channelNumber = other.channelNumber; + } + return *this; + } +}; + +class GPRDataModel { +public: + struct Header { + int numTraces = 0; // LAST TRACE + int samplesPerTrace = 0; // SAMPLES + double timeWindowNs = 0.0; // TIMEWINDOW (ns) + double distanceInc = 0.0; // DISTANCE INTERVAL (m) + double antennaFreq = 0.0; // FREQUENCY (MHz) + QString antennaType; // ANTENNAS + QString date; // DATE + QString timeStr; // TIME + // 三维相关参数 + int numberOfChannels = 1; // NUMBER_OF_CH + QVector chXOffsets; // CH_X_OFFSETS + QVector chYOffsets; // CH_Y_OFFSETS + QString units; // UNITS + double startPosition = 0.0; // START POSITION + double stopPosition = 0.0; // STOP POSITION + double waveVelocity = 0.1; // 波速，用于深度计算 + double timeIntervalNs = 0.0; // TIME INTERVAL (ns)，用于计算时窗 + + // 原始映射，用于调试或扩展 + QMap rawParams; + }; + + Header header; + QVector traces; + + // 三维数据相关属性 + int tracesPerChannel = 0; // 每个通道的道数 + int channels = 0; // 实际通道数 + double totalDistance = 0.0; // 总距离 + + // 三维数据体 - 存储为[channel][trace][sample]格式 + QVector>> volumeData; // [channel][trace][sample] + + bool isEmpty() const { return traces.isEmpty(); } + void clear() { + traces.clear(); + header = Header{}; + tracesPerChannel = 0; + channels = 0; + totalDistance = 0.0; + volumeData.clear(); + } + + // 获取指定通道的道数 + int getTracesPerChannel() const { + if (header.numTraces <= 0) { + return 0; + } + if (header.numberOfChannels > 0) { + return header.numTraces / header.numberOfChannels; + } + return header.numTraces; + } + + int getTraceIndex(int channel, int traceInChannel) const { + const int nChannels = header.numberOfChannels > 0 ? header.numberOfChannels : 1; + if (channel < 0 || channel >= nChannels || + traceInChannel < 0 || traceInChannel >= getTracesPerChannel()) { + return -1; + } + return traceInChannel * nChannels + channel; + } + + // 获取指定通道和道号的数据 + const RadarTrace* getTrace(int channel, int traceInChannel) const { + const int index = getTraceIndex(channel, traceInChannel); + if (index >= 0 && index < traces.size()) { + return &traces[index]; + } + return nullptr; + } + + // 计算深度（米） + double calculateDepth(int sampleIndex) const { + if (header.samplesPerTrace > 0) { + // 使用默认波速0.1m/ns，除以2是因为往返时间 + double timePerSample = header.timeWindowNs / header.samplesPerTrace; + double timeNs = sampleIndex * timePerSample; + return timeNs * header.waveVelocity / 2.0; // 深度 = 时间 * 波速 / 2 + } + return 0.0; + } + + // 计算距离（米） + double calculateDistance(int traceIndex) const { + if (header.distanceInc > 0) { + return traceIndex * header.distanceInc; + } + return 0.0; + } + + // 构建三维数据体 + void buildVolumeData() { + if (traces.isEmpty()) return; + + int nChannels = header.numberOfChannels > 0 ? header.numberOfChannels : 1; + int nTracesPerChannel = getTracesPerChannel(); + int nSamples = header.samplesPerTrace; + + volumeData.resize(nChannels); + for (int c = 0; c < nChannels; ++c) { + volumeData[c].resize(nTracesPerChannel); + for (int t = 0; t < nTracesPerChannel; ++t) { + volumeData[c][t].resize(nSamples); + } + } + + // 填充数据 + for (int i = 0; i < traces.size(); ++i) { + int channel = i % nChannels; + int traceInChannel = i / nChannels; + + if (channel < volumeData.size() && traceInChannel < volumeData[channel].size()) { + for (int s = 0; s < traces[i].amplitudes.size() && s < nSamples; ++s) { + volumeData[channel][traceInChannel][s] = traces[i].amplitudes[s]; + } + } + } + } + + // 获取三维数据点 + short getVolumeValue(int channel, int trace, int sample) const { + if (channel < 0 || channel >= volumeData.size() || + trace < 0 || trace >= volumeData[channel].size() || + sample < 0 || sample >= volumeData[channel][trace].size()) { + return 0; + } + return volumeData[channel][trace][sample]; + } + + // 获取全局最小最大值 + void getGlobalMinMax(short& minVal, short& maxVal) const { + if (volumeData.isEmpty()) { + minVal = 0; + maxVal = 0; + return; + } + + minVal = 32767; + maxVal = -32768; + + for (const auto& channelData : volumeData) { + for (const auto& traceData : channelData) { + for (short val : traceData) { + if (val < minVal) minVal = val; + if (val > maxVal) maxVal = val; + } + } + } + } + + // 获取通道数 + int getChannelCount() const { + if (!volumeData.isEmpty()) return volumeData.size(); + return qMax(1, header.numberOfChannels); + } + + // 获取每通道的迹数 + int getTraceCountPerChannel() const { + if (!volumeData.isEmpty() && !volumeData[0].isEmpty()) return volumeData[0].size(); + int nCh = qMax(1, header.numberOfChannels); + return header.numTraces / nCh; + } + + // 获取每迹的样本数 + int getSampleCount() const { + if (!volumeData.isEmpty() && !volumeData[0].isEmpty()) return volumeData[0][0].size(); + return header.samplesPerTrace; + } + + GPRDataModel() = default; + GPRDataModel(const GPRDataModel&) = default; + GPRDataModel& operator=(const GPRDataModel&) = default; + GPRDataModel(GPRDataModel&&) = default; + GPRDataModel& operator=(GPRDataModel&&) = default; +}; + +struct TrajectoryEditPoint { + QVector3D rawPosition; + QVector3D editedPosition; + bool enabled = true; + bool manuallyEdited = false; + bool distanceOutlier = false; + bool speedOutlier = false; + bool angleOutlier = false; +}; + +struct TrajectoryEditSettings { + bool distanceFilterEnabled = true; + double maxSegmentDistanceM = 1.0; + bool speedFilterEnabled = false; + double maxSpeedMps = 5.0; + double traceIntervalSec = 0.05; + bool angleFilterEnabled = true; + double maxTurnAngleDeg = 60.0; + int interpolationMode = 0; // 0 linear, 1 spline + bool preserveManualEdits = true; +}; + +// 单条测线结构（含GPS） +struct SurveyLine { + QString name; // 测线编号，如 "_000", "_001" + QString displayName; // 显示名称，如 "碧新路_000" + QString radFilePath; // .rad 文件完整路径 + GPRDataModel data; // 雷达数据 + QVector gpsPositions; // 每个trace的GPS坐标 (X,Y,Z) + QVector trajectoryEditPoints; // 可编辑中心线轨迹点 + TrajectoryEditSettings trajectoryEditSettings; + + // 三维轨迹相关（新增） + SurveyGeometry geometry; // 天线几何参数 + QVector> channelTrajectories; // [channel][trace] 绝对坐标 + + bool hasGps() const { return !gpsPositions.isEmpty(); } + bool hasValidTrajectories() const { return !channelTrajectories.isEmpty() && !channelTrajectories[0].isEmpty(); } + bool isEmpty() const { return data.isEmpty(); } + void clear() { + name.clear(); + displayName.clear(); + radFilePath.clear(); + data.clear(); + gpsPositions.clear(); + trajectoryEditPoints.clear(); + trajectoryEditSettings = TrajectoryEditSettings{}; + geometry = SurveyGeometry{}; + channelTrajectories.clear(); + } + + SurveyLine() = default; + SurveyLine(const SurveyLine&) = default; + SurveyLine& operator=(const SurveyLine&) = default; + SurveyLine(SurveyLine&&) = default; + SurveyLine& operator=(SurveyLine&&) = default; +}; + +Q_DECLARE_METATYPE(SurveyLine) + +#endif // GPRDATAMODEL_H \ No newline at end of file diff --git a/external/gpr3dviewer/ImpulseMultiChannelConverter.cpp b/external/gpr3dviewer/ImpulseMultiChannelConverter.cpp new file mode 100644 index 0000000..1527992 --- /dev/null +++ b/external/gpr3dviewer/ImpulseMultiChannelConverter.cpp @@ -0,0 +1,406 @@ +#include "ImpulseMultiChannelConverter.h" + +#include "IprhParser.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +bool ImpulseMultiChannelConverter::isMultiChannelImpulseHeader(const QString &headerPath, + QString *dirPath, + QString *baseName) +{ + QFileInfo fi(headerPath); + QRegularExpression reMultiChannel(QStringLiteral("^(.*)_A(\\d+)$"), + QRegularExpression::CaseInsensitiveOption); + QRegularExpressionMatch match = reMultiChannel.match(fi.completeBaseName()); + if (!match.hasMatch()) + return false; + + const QString surveyBase = match.captured(1); + QDir dir(fi.absolutePath()); + const QStringList channels = dir.entryList(QStringList() << surveyBase + QStringLiteral("_A*.iprh"), + QDir::Files); + if (channels.size() <= 1) + return false; + + if (dirPath) + *dirPath = fi.absolutePath(); + if (baseName) + *baseName = surveyBase; + return true; +} + +bool ImpulseMultiChannelConverter::buildPlan(const QString &dirPath, + const QString &baseName, + ConversionPlan &plan, + QString *errorMessage) +{ + plan = ConversionPlan{}; + plan.dirPath = dirPath; + plan.baseName = baseName; + + QDir dir(dirPath); + dir.setFilter(QDir::Files | QDir::NoDotAndDotDot); + + QRegularExpression reChannel(QStringLiteral("^%1_A(\\d+)\\.iprh$").arg(QRegularExpression::escape(baseName)), + QRegularExpression::CaseInsensitiveOption); + QMap channelHeaders; + for (const QFileInfo &fi : dir.entryInfoList()) { + QRegularExpressionMatch match = reChannel.match(fi.fileName()); + if (!match.hasMatch()) + continue; + + const int chNum = match.captured(1).toInt(); + const QString iprhPath = fi.absoluteFilePath(); + QString iprbPath = iprhPath; + iprbPath.replace(QStringLiteral(".iprh"), QStringLiteral(".iprb"), Qt::CaseInsensitive); + if (QFile::exists(iprbPath)) { + channelHeaders.insert(chNum, iprhPath); + } else { + qDebug() << "Missing .iprb for" << iprhPath; + } + } + + if (channelHeaders.isEmpty()) { + if (errorMessage) + *errorMessage = QStringLiteral("未找到 Impulse 多通道文件：%1").arg(baseName); + return false; + } + + GPRDataModel::Header masterHeader; + if (!IprhParser::parseHeaderOnly(channelHeaders.first(), masterHeader)) { + if (errorMessage) + *errorMessage = QStringLiteral("无法解析主通道头文件：%1").arg(channelHeaders.first()); + return false; + } + + QVector xOffsets; + QVector yOffsets; + xOffsets.reserve(channelHeaders.size()); + yOffsets.reserve(channelHeaders.size()); + + qint64 minTraceCount = std::numeric_limits::max(); + double antennaSeparation = 0.0; + + for (auto it = channelHeaders.begin(); it != channelHeaders.end(); ++it) { + ChannelInfo info; + info.channelNumber = it.key(); + info.iprhPath = it.value(); + info.iprbPath = info.iprhPath; + info.iprbPath.replace(QStringLiteral(".iprh"), QStringLiteral(".iprb"), Qt::CaseInsensitive); + + if (!IprhParser::parseHeaderOnly(info.iprhPath, info.header)) { + if (errorMessage) + *errorMessage = QStringLiteral("无法解析通道头文件：%1").arg(info.iprhPath); + return false; + } + + if (info.header.samplesPerTrace != masterHeader.samplesPerTrace) { + if (errorMessage) + *errorMessage = QStringLiteral("多通道 SAMPLES 不一致：%1").arg(info.iprhPath); + return false; + } + if (!qFuzzyCompare(info.header.timeIntervalNs, masterHeader.timeIntervalNs) && + info.header.timeIntervalNs > 0 && masterHeader.timeIntervalNs > 0) { + qDebug() << "Warning: Inconsistent TIME INTERVAL across channels" << info.iprhPath; + } + if (!qFuzzyCompare(info.header.distanceInc, masterHeader.distanceInc) && + info.header.distanceInc > 0 && masterHeader.distanceInc > 0) { + qDebug() << "Warning: Inconsistent DISTANCE INTERVAL across channels" << info.iprhPath; + } + + if (info.header.rawParams.contains(QStringLiteral("CH_X_OFFSET"))) { + info.xOffset = info.header.rawParams.value(QStringLiteral("CH_X_OFFSET")).toFloat(); + } else if (info.header.rawParams.contains(QStringLiteral("CH_OFFSET_X"))) { + info.xOffset = info.header.rawParams.value(QStringLiteral("CH_OFFSET_X")).toFloat(); + } else if (!info.header.chXOffsets.isEmpty()) { + info.xOffset = info.header.chXOffsets.first(); + } + if (info.header.rawParams.contains(QStringLiteral("CH_Y_OFFSET"))) { + info.yOffset = info.header.rawParams.value(QStringLiteral("CH_Y_OFFSET")).toFloat(); + } else if (info.header.rawParams.contains(QStringLiteral("CH_OFFSET_Y"))) { + info.yOffset = info.header.rawParams.value(QStringLiteral("CH_OFFSET_Y")).toFloat(); + } else if (!info.header.chYOffsets.isEmpty()) { + info.yOffset = info.header.chYOffsets.first(); + } + xOffsets.append(info.xOffset); + yOffsets.append(info.yOffset); + + if (info.header.rawParams.contains(QStringLiteral("ANTENNA SEPARATION"))) { + antennaSeparation = info.header.rawParams.value(QStringLiteral("ANTENNA SEPARATION")).toDouble(); + } + + const qint64 traceBytes = static_cast(info.header.samplesPerTrace) * sizeof(qint16); + if (traceBytes <= 0) { + if (errorMessage) + *errorMessage = QStringLiteral("通道采样数无效：%1").arg(info.iprhPath); + return false; + } + info.traceCount = QFileInfo(info.iprbPath).size() / traceBytes; + if (info.traceCount <= 0) { + if (errorMessage) + *errorMessage = QStringLiteral("通道无有效数据：%1").arg(info.iprbPath); + return false; + } + minTraceCount = std::min(minTraceCount, info.traceCount); + plan.channels.append(info); + } + + if (plan.channels.isEmpty() || minTraceCount <= 0 || minTraceCount == std::numeric_limits::max()) { + if (errorMessage) + *errorMessage = QStringLiteral("Impulse 多通道没有有效 trace：%1").arg(baseName); + return false; + } + + for (const ChannelInfo &info : plan.channels) { + if (info.traceCount != minTraceCount) { + qDebug() << "Warning: Inconsistent trace count across channels. Using minimum:" + << minTraceCount << "channel has:" << info.traceCount << info.iprhPath; + } + } + + bool allZeroOffsets = true; + for (float value : xOffsets) { + if (!qFuzzyIsNull(value)) { + allZeroOffsets = false; + break; + } + } + if (allZeroOffsets && antennaSeparation > 1e-6 && plan.channels.size() > 1) { + const double totalWidth = (plan.channels.size() - 1) * antennaSeparation; + const double startX = -totalWidth / 2.0; + for (int i = 0; i < xOffsets.size(); ++i) { + xOffsets[i] = static_cast(startX + i * antennaSeparation); + plan.channels[i].xOffset = xOffsets[i]; + } + qDebug() << "Computed symmetric X offsets from antenna separation:" << antennaSeparation; + } + + plan.channelCount = plan.channels.size(); + plan.tracesPerChannel = static_cast(minTraceCount); + plan.samplesPerTrace = masterHeader.samplesPerTrace; + plan.traceByteSize = plan.samplesPerTrace * static_cast(sizeof(qint16)); + plan.totalOutputTraces = static_cast(plan.tracesPerChannel) * plan.channelCount; + plan.outputHeader = masterHeader; + plan.outputHeader.numberOfChannels = plan.channelCount; + plan.outputHeader.numTraces = static_cast(plan.totalOutputTraces); + plan.outputHeader.chXOffsets = xOffsets; + plan.outputHeader.chYOffsets = yOffsets; + if (plan.outputHeader.timeWindowNs <= 0.0 && plan.outputHeader.timeIntervalNs > 0.0) { + plan.outputHeader.timeWindowNs = plan.outputHeader.samplesPerTrace * plan.outputHeader.timeIntervalNs; + } + + const QString basePath = dir.absoluteFilePath(baseName + QStringLiteral("_mala_converted")); + plan.outputRadPath = basePath + QStringLiteral(".rad"); + plan.outputRd3Path = basePath + QStringLiteral(".rd3"); + return true; +} + +bool ImpulseMultiChannelConverter::convertStreaming(const ConversionPlan &plan, + const Options &options, + QString *radFilePath, + QString *errorMessage, + CancelFn cancel, + ProgressFn progress) +{ + if (plan.channelCount <= 0 || plan.tracesPerChannel <= 0 || plan.traceByteSize <= 0) { + if (errorMessage) + *errorMessage = QStringLiteral("Impulse 转换计划无效"); + return false; + } + + if (options.reuseExistingIfValid && convertedFilesAreValid(plan)) { + if (radFilePath) + *radFilePath = plan.outputRadPath; + return true; + } + + if (!options.overwriteExisting && (QFile::exists(plan.outputRadPath) || QFile::exists(plan.outputRd3Path))) { + if (errorMessage) + *errorMessage = QStringLiteral("转换输出文件已存在：%1").arg(plan.outputRadPath); + return false; + } + + const QString tmpRd3Path = plan.outputRd3Path + QStringLiteral(".tmp"); + const QString tmpRadPath = plan.outputRadPath + QStringLiteral(".tmp"); + QFile::remove(tmpRd3Path); + QFile::remove(tmpRadPath); + + QVector> inputFiles; + inputFiles.reserve(plan.channels.size()); + for (const ChannelInfo &channel : plan.channels) { + auto file = QSharedPointer::create(channel.iprbPath); + if (!file->open(QIODevice::ReadOnly)) { + if (errorMessage) + *errorMessage = QStringLiteral("无法读取 Impulse 通道数据：%1").arg(channel.iprbPath); + return false; + } + inputFiles.append(file); + } + + QFile rd3File(tmpRd3Path); + if (!rd3File.open(QIODevice::WriteOnly | QIODevice::Truncate)) { + if (errorMessage) + *errorMessage = QStringLiteral("无法写入转换后的 RD3 文件：%1").arg(tmpRd3Path); + return false; + } + + const qint64 bytesPerTracePosition = plan.traceByteSize * plan.channelCount; + const qint64 chunkBudget = qMax(bytesPerTracePosition, options.maxChunkBytes); + const int chunkTraces = qMax(1, chunkBudget / bytesPerTracePosition); + QVector channelBuffers(plan.channelCount); + + qint64 tracesDone = 0; + while (tracesDone < plan.tracesPerChannel) { + if (cancel && cancel()) { + rd3File.close(); + QFile::remove(tmpRd3Path); + QFile::remove(tmpRadPath); + if (errorMessage) + *errorMessage = QStringLiteral("Impulse 多通道转换已取消"); + return false; + } + + const int currentChunk = qMin(chunkTraces, plan.tracesPerChannel - tracesDone); + const qint64 expectedChannelBytes = currentChunk * plan.traceByteSize; + for (int ch = 0; ch < plan.channelCount; ++ch) { + channelBuffers[ch] = inputFiles[ch]->read(expectedChannelBytes); + if (channelBuffers[ch].size() != expectedChannelBytes) { + rd3File.close(); + QFile::remove(tmpRd3Path); + if (errorMessage) { + *errorMessage = QStringLiteral("读取通道数据不完整：%1").arg(plan.channels[ch].iprbPath); + } + return false; + } + } + + for (int localTrace = 0; localTrace < currentChunk; ++localTrace) { + const qint64 offset = localTrace * plan.traceByteSize; + for (int ch = 0; ch < plan.channelCount; ++ch) { + const qint64 written = rd3File.write(channelBuffers[ch].constData() + offset, plan.traceByteSize); + if (written != plan.traceByteSize) { + rd3File.close(); + QFile::remove(tmpRd3Path); + if (errorMessage) + *errorMessage = QStringLiteral("写入转换后的 RD3 文件失败：%1").arg(tmpRd3Path); + return false; + } + } + } + + tracesDone += currentChunk; + if (progress) { + progress(Progress{tracesDone, plan.tracesPerChannel, + QStringLiteral("正在转换 Impulse 多通道 %1/%2") + .arg(tracesDone) + .arg(plan.tracesPerChannel)}); + } + } + + rd3File.close(); + if (rd3File.error() != QFile::NoError) { + QFile::remove(tmpRd3Path); + if (errorMessage) + *errorMessage = QStringLiteral("保存转换后的 RD3 文件失败：%1").arg(rd3File.errorString()); + return false; + } + + ConversionPlan tmpPlan = plan; + tmpPlan.outputRadPath = tmpRadPath; + if (!writeRadHeader(tmpRadPath, tmpPlan, errorMessage)) { + QFile::remove(tmpRd3Path); + QFile::remove(tmpRadPath); + return false; + } + + QFile::remove(plan.outputRd3Path); + QFile::remove(plan.outputRadPath); + if (!QFile::rename(tmpRd3Path, plan.outputRd3Path)) { + QFile::remove(tmpRd3Path); + QFile::remove(tmpRadPath); + if (errorMessage) + *errorMessage = QStringLiteral("无法替换转换后的 RD3 文件：%1").arg(plan.outputRd3Path); + return false; + } + if (!QFile::rename(tmpRadPath, plan.outputRadPath)) { + QFile::remove(plan.outputRd3Path); + QFile::remove(tmpRadPath); + if (errorMessage) + *errorMessage = QStringLiteral("无法替换转换后的 RAD 文件：%1").arg(plan.outputRadPath); + return false; + } + + if (radFilePath) + *radFilePath = plan.outputRadPath; + + qDebug() << "Impulse multi-channel streaming conversion OK:" << plan.outputRadPath << plan.outputRd3Path; + return true; +} + +bool ImpulseMultiChannelConverter::writeRadHeader(const QString &radFilePath, + const ConversionPlan &plan, + QString *errorMessage) +{ + QFile radFile(radFilePath); + if (!radFile.open(QIODevice::WriteOnly | QIODevice::Text | QIODevice::Truncate)) { + if (errorMessage) + *errorMessage = QStringLiteral("无法写入转换后的 RAD 文件：%1").arg(radFilePath); + return false; + } + + QTextStream rad(&radFile); + const auto &header = plan.outputHeader; + rad << "# Converted from Impulse multi-channel survey: " << plan.baseName << '\n'; + rad << "SAMPLES: " << header.samplesPerTrace << '\n'; + rad << "LAST TRACE: " << plan.totalOutputTraces << '\n'; + rad << "TIMEWINDOW: " << QLocale::c().toString(header.timeWindowNs, 'g', 12) << '\n'; + rad << "TIME INTERVAL: " << QLocale::c().toString(header.timeIntervalNs, 'g', 12) << '\n'; + rad << "DISTANCE INTERVAL: " << QLocale::c().toString(header.distanceInc, 'g', 12) << '\n'; + rad << "ANTENNAS: " << QLocale::c().toString(header.antennaFreq, 'g', 12) << '\n'; + if (!header.antennaType.isEmpty()) rad << "ANTENNA: " << header.antennaType << '\n'; + if (!header.date.isEmpty()) rad << "DATE: " << header.date << '\n'; + if (!header.timeStr.isEmpty()) rad << "TIME: " << header.timeStr << '\n'; + rad << "NUMBER_OF_CH: " << qMax(1, header.numberOfChannels) << '\n'; + rad << "CH_X_OFFSETS: " << formatFloatList(header.chXOffsets) << '\n'; + rad << "CH_Y_OFFSETS: " << formatFloatList(header.chYOffsets) << '\n'; + if (!header.units.isEmpty()) rad << "UNITS: " << header.units << '\n'; + rad << "START POSITION: " << QLocale::c().toString(header.startPosition, 'g', 12) << '\n'; + const double stopPosition = header.stopPosition > header.startPosition + ? header.stopPosition + : header.startPosition + plan.tracesPerChannel * header.distanceInc; + rad << "STOP POSITION: " << QLocale::c().toString(stopPosition, 'g', 12) << '\n'; + radFile.close(); + return true; +} + +bool ImpulseMultiChannelConverter::convertedFilesAreValid(const ConversionPlan &plan) +{ + QFileInfo radInfo(plan.outputRadPath); + QFileInfo rd3Info(plan.outputRd3Path); + if (!radInfo.exists() || !rd3Info.exists()) + return false; + + const qint64 expectedBytes = plan.totalOutputTraces * plan.traceByteSize; + return rd3Info.size() == expectedBytes; +} + +QString ImpulseMultiChannelConverter::formatFloatList(const QVector &values) +{ + QStringList parts; + parts.reserve(values.size()); + for (float value : values) { + parts.append(QLocale::c().toString(value, 'g', 10)); + } + return parts.join(' '); +} diff --git a/external/gpr3dviewer/ImpulseMultiChannelConverter.h b/external/gpr3dviewer/ImpulseMultiChannelConverter.h new file mode 100644 index 0000000..dde2fa9 --- /dev/null +++ b/external/gpr3dviewer/ImpulseMultiChannelConverter.h @@ -0,0 +1,75 @@ +#ifndef IMPULSEMULTICHANNELCONVERTER_H +#define IMPULSEMULTICHANNELCONVERTER_H + +#include "GPRDataModel.h" + +#include +#include +#include + +class ImpulseMultiChannelConverter { +public: + struct ChannelInfo { + int channelNumber = 0; + QString iprhPath; + QString iprbPath; + GPRDataModel::Header header; + qint64 traceCount = 0; + float xOffset = 0.0f; + float yOffset = 0.0f; + }; + + struct ConversionPlan { + QString dirPath; + QString baseName; + QVector channels; + GPRDataModel::Header outputHeader; + int channelCount = 0; + int tracesPerChannel = 0; + qint64 totalOutputTraces = 0; + qint64 samplesPerTrace = 0; + qint64 traceByteSize = 0; + QString outputRadPath; + QString outputRd3Path; + }; + + struct Options { + qint64 maxChunkBytes = 32 * 1024 * 1024; + bool overwriteExisting = true; + bool reuseExistingIfValid = false; + }; + + struct Progress { + qint64 tracesDone = 0; + qint64 tracesTotal = 0; + QString message; + }; + + using CancelFn = std::function; + using ProgressFn = std::function; + + static bool isMultiChannelImpulseHeader(const QString &headerPath, + QString *dirPath = nullptr, + QString *baseName = nullptr); + + static bool buildPlan(const QString &dirPath, + const QString &baseName, + ConversionPlan &plan, + QString *errorMessage = nullptr); + + static bool convertStreaming(const ConversionPlan &plan, + const Options &options, + QString *radFilePath = nullptr, + QString *errorMessage = nullptr, + CancelFn cancel = {}, + ProgressFn progress = {}); + +private: + static bool writeRadHeader(const QString &radFilePath, + const ConversionPlan &plan, + QString *errorMessage); + static bool convertedFilesAreValid(const ConversionPlan &plan); + static QString formatFloatList(const QVector &values); +}; + +#endif // IMPULSEMULTICHANNELCONVERTER_H diff --git a/external/gpr3dviewer/IprhParser.cpp b/external/gpr3dviewer/IprhParser.cpp new file mode 100644 index 0000000..608f6ca --- /dev/null +++ b/external/gpr3dviewer/IprhParser.cpp @@ -0,0 +1,625 @@ +#include "IprhParser.h" +#include "PerformanceLogger.h" + +#include "ImpulseMultiChannelConverter.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/** + * @brief 加载 Impulse 系列雷达 IPRH 数据（iprh头文件 + iprb纯二进制波形） + * @param iprhFilePath .iprh 文本配置头文件路径 + * @param model 输出 GPR 全局数据模型 + * @return true 加载解析成功；false 失败 + * @note 存储结构：.iprh 存储全部仪器/测线参数；.iprb 无文件头，从头到尾全是 short16 振幅采样 + */ +bool IprhParser::loadFromIprh(const QString &iprhFilePath, GPRDataModel &model) { + SCOPED_PERF_TIMER("Parser", "IprhParser::loadFromIprh"); + + model.clear(); + model.header = GPRDataModel::Header{}; + + // 第一步：解析 iprh 文本头 + if (!parseIprhHeader(iprhFilePath, model.header)) { + qDebug() << "Error: Failed to parse .iprh header file:" << iprhFilePath; + return false; + } + + // 自动匹配同目录同名二进制数据文件 .iprb + QFileInfo iprhInfo(iprhFilePath); + QString iprbPath = iprhInfo.absolutePath() + "/" + iprhInfo.completeBaseName() + ".iprb"; + + if (!QFile::exists(iprbPath)) { + qDebug() << "Error: Matching .iprb binary file not found at:" << iprbPath; + return false; + } + + // 如果头文件没有 LAST TRACE，从二进制文件大小推算 + if (model.header.numTraces <= 0) { + QFile binaryFile(iprbPath); + if (binaryFile.open(QIODevice::ReadOnly)) { + qint64 fileSize = binaryFile.size(); + qint64 traceBytes = static_cast(model.header.samplesPerTrace) * sizeof(short); + if (traceBytes > 0) { + model.header.numTraces = static_cast(fileSize / traceBytes); + qDebug() << "Inferred numTraces from binary size:" << model.header.numTraces; + } + binaryFile.close(); + } + } + if (model.header.timeWindowNs <= 0.0 && model.header.timeIntervalNs > 0.0) { + model.header.timeWindowNs = model.header.samplesPerTrace * model.header.timeIntervalNs; + qDebug() << "Inferred timeWindowNs from samples * timeInterval:" << model.header.timeWindowNs; + } + + if (model.header.numTraces <= 0) { + qDebug() << "Error: Unable to determine numTraces from header or binary file"; + return false; + } + + // 第二步：读取纯二进制波形 + return loadIprbBinary(iprbPath, model); +} + +/** + * @brief 解析 IPRH 文本头文件，提取雷达采集关键参数 + * @param iprhFilePath iprh 文本文件路径 + * @param header 待填充头部参数结构体 + * @return 解析成功返回 true + */ +bool IprhParser::parseHeaderOnly(const QString &iprhFilePath, GPRDataModel::Header &header) +{ + return parseIprhHeader(iprhFilePath, header); +} + +bool IprhParser::parseIprhHeader(const QString &iprhFilePath, GPRDataModel::Header &header) { + SCOPED_PERF_TIMER("Parser", "IprhParser::parseIprhHeader"); + + QFile file(iprhFilePath); + if (!file.open(QIODevice::ReadOnly | QIODevice::Text)) { + qDebug() << "Error: Cannot open iprh file for read"; + return false; + } + + QTextStream in(&file); + + while (!in.atEnd()) { + QString line = in.readLine().trimmed(); + if (line.isEmpty()) continue; + + int sepPos = line.indexOf(':'); + if (sepPos == -1) continue; + + QString key = line.left(sepPos).trimmed(); + QString value = line.mid(sepPos + 1).trimmed(); + + header.rawParams[key] = value; + + if (key == "SAMPLES") { + header.samplesPerTrace = extractInt(value); + } else if (key == "LAST TRACE") { + header.numTraces = extractInt(value); + } else if (key == "TIMEWINDOW") { + header.timeWindowNs = extractDouble(value); + } else if (key == "DISTANCE INTERVAL") { + header.distanceInc = extractDouble(value); + } else if (key == "TIME INTERVAL") { + header.timeIntervalNs = extractDouble(value); + } else if (key == "ANTENNAS") { + header.antennaFreq = extractDouble(value); + } else if (key == "ANTENNA") { + header.antennaType = value; + } else if (key == "DATE") { + header.date = value; + } else if (key == "START TIME" || key == "TIME") { + header.timeStr = value; + } else if (key == "CHANNELS" || key == "NUMBER_OF_CH") { + header.numberOfChannels = extractInt(value); + } else if (key == "CH_X_OFFSET" || key == "CH_OFFSET_X") { + // 单通道偏移量（Impulse 单通道文件） + if (!value.isEmpty()) { + header.chXOffsets.append(value.toFloat()); + } + } else if (key == "CH_Y_OFFSET" || key == "CH_OFFSET_Y") { + if (!value.isEmpty()) { + header.chYOffsets.append(value.toFloat()); + } + } else if (key == "CH_X_OFFSETS") { + // 兼容 Mala Mira 风格多值空格分隔 + QStringList offsets = value.split(' ', Qt::SkipEmptyParts); + for (const QString &offset : offsets) { + header.chXOffsets.append(offset.toFloat()); + } + } else if (key == "CH_Y_OFFSETS") { + QStringList offsets = value.split(' ', Qt::SkipEmptyParts); + for (const QString &offset : offsets) { + header.chYOffsets.append(offset.toFloat()); + } + } else if (key == "UNITS") { + header.units = value; + } else if (key == "START POSITION") { + header.startPosition = extractDouble(value); + } else if (key == "STOP POSITION") { + header.stopPosition = extractDouble(value); + } + } + + file.close(); + + if (header.samplesPerTrace <= 0) { + qDebug() << "Error: Invalid SAMPLES value in iprh file"; + return false; + } + + header.waveVelocity = 0.1; + + qDebug() << "==== IPRH Header Parse Complete ====" + << "\nSamples per trace:" << header.samplesPerTrace + << "\nTotal traces:" << header.numTraces + << "\nTime window(ns):" << header.timeWindowNs + << "\nChannel count:" << header.numberOfChannels + << "\nX offset array size:" << header.chXOffsets.size() + << "\nY offset array size:" << header.chYOffsets.size(); + + return true; +} + +/** + * @brief 读取纯 IPRB 二进制数据（无任何文件头，全连续 short16 振幅） + * @param iprbFilePath iprb 波形文件路径 + * @param model 绑定头部参数，填充完整 traces 波形数组 + * @return 读取成功 true + */ +bool IprhParser::loadIprbBinary(const QString &iprbFilePath, GPRDataModel &model) { + SCOPED_PERF_TIMER("Parser", "IprhParser::loadIprbBinary"); + + QFile file(iprbFilePath); + if (!file.open(QIODevice::ReadOnly)) { + qDebug() << "Error: Open iprb binary failed:" << iprbFilePath; + return false; + } + + const int samplesPerTrace = model.header.samplesPerTrace; + const int totalTraceCount = model.header.numTraces; + + const qint64 dataStartOffset = 0; + file.seek(dataStartOffset); + + const qint64 singleTraceByteSize = samplesPerTrace * sizeof(short); + const qint64 fullExpectedBytes = totalTraceCount * singleTraceByteSize; + const qint64 realFileBytes = file.size(); + + if (realFileBytes < fullExpectedBytes) { + qDebug() << "Warning: IPRB file size smaller than theoretical data size!" + << "Expected:" << fullExpectedBytes << "Actual:" << realFileBytes; + } + + model.traces.reserve(totalTraceCount); + + QByteArray traceBuffer; + traceBuffer.resize(singleTraceByteSize); + + const int channelCnt = model.header.numberOfChannels > 0 ? model.header.numberOfChannels : 1; + model.channels = channelCnt; + model.tracesPerChannel = totalTraceCount / channelCnt; + + if (model.header.distanceInc > 1e-6) { + model.totalDistance = static_cast(model.tracesPerChannel * model.header.distanceInc); + } else { + model.totalDistance = static_cast(model.header.stopPosition - model.header.startPosition); + } + + for (int traceGlobalIdx = 0; traceGlobalIdx < totalTraceCount; ++traceGlobalIdx) { + if (file.atEnd()) { + qDebug() << "Warning: File ended early at global trace index" << traceGlobalIdx; + break; + } + + qint64 readBytes = file.read(traceBuffer.data(), traceBuffer.size()); + if (readBytes != traceBuffer.size()) { + qDebug() << "Warning: Trace" << traceGlobalIdx << "incomplete byte read"; + break; + } + + RadarTrace oneTrace; + oneTrace.amplitudes.resize(samplesPerTrace); + const short* rawShortBuf = reinterpret_cast(traceBuffer.constData()); + + for (int s = 0; s < samplesPerTrace; s++) { + oneTrace.amplitudes[s] = rawShortBuf[s]; + } + + int chNo = traceGlobalIdx % channelCnt; + int traceInChIdx = traceGlobalIdx / channelCnt; + oneTrace.channelNumber = chNo; + + float xOff = 0.0f; + float yOff = 0.0f; + if (chNo < model.header.chXOffsets.size()) xOff = model.header.chXOffsets[chNo]; + if (chNo < model.header.chYOffsets.size()) yOff = model.header.chYOffsets[chNo]; + + float lineDist = static_cast(model.header.startPosition + traceInChIdx * model.header.distanceInc); + oneTrace.position = QVector3D(xOff, lineDist - yOff, 0.0f); + + model.traces.append(std::move(oneTrace)); + } + + file.close(); + + if (model.traces.isEmpty()) { + qDebug() << "Error: No valid traces loaded from iprb"; + return false; + } + + qDebug() << "IPRB Binary Load OK, total valid traces:" << model.traces.size(); + return true; +} + +/** + * @brief 加载 Impulse 多通道数据（每个通道一个 .iprh + .iprb） + * @param dirPath 数据所在目录 + * @param baseName 测线基础名（如 "明星路_001"） + * @param model 输出合并后的 GPR 数据模型 + * @return true 加载合并成功 + * + * 文件命名约定：baseName_A01.iprh / .iprb ... baseName_A14.iprh / .iprb + * 合并后 traces 按 Mala Mira 格式交错：trace0=ch0_pos0, trace1=ch1_pos0, ... + */ +bool IprhParser::convertImpulseMultiChannelToMala(const QString &dirPath, + const QString &baseName, + QString *radFilePath, + QString *errorMessage) +{ + ImpulseMultiChannelConverter::ConversionPlan plan; + if (!ImpulseMultiChannelConverter::buildPlan(dirPath, baseName, plan, errorMessage)) { + if (errorMessage && errorMessage->isEmpty()) { + *errorMessage = QStringLiteral("多通道 Impulse 数据合并失败：%1").arg(baseName); + } + return false; + } + + ImpulseMultiChannelConverter::Options options; + options.overwriteExisting = true; + options.reuseExistingIfValid = false; + if (!ImpulseMultiChannelConverter::convertStreaming(plan, options, radFilePath, errorMessage)) { + return false; + } + + qDebug() << "Impulse multi-channel converted to Mala Mira files:" << plan.outputRadPath << plan.outputRd3Path; + return true; +} + +/** + * @brief 加载 Impulse 多通道数据（每个通道一个 .iprh + .iprb） + * @param dirPath 数据所在目录 + * @param baseName 测线基础名（如 "明星路_001"） + * @param model 输出合并后的 GPR 数据模型 + * @return true 加载合并成功 + * + * 文件命名约定：baseName_A01.iprh / .iprb ... baseName_A14.iprh / .iprb + * 合并后 traces 按 Mala Mira 格式交错：trace0=ch0_pos0, trace1=ch1_pos0, ... + */ +bool IprhParser::loadImpulseMultiChannel(const QString &dirPath, + const QString &baseName, + GPRDataModel &model) +{ + SCOPED_PERF_TIMER("Parser", "IprhParser::loadImpulseMultiChannel"); + + model.clear(); + model.header = GPRDataModel::Header{}; + + // 1. 发现所有通道文件并排序 + QDir dir(dirPath); + dir.setFilter(QDir::Files | QDir::NoDotAndDotDot); + + QRegularExpression reChannel(QStringLiteral("^%1_A(\\d+)\\.iprh$").arg(QRegularExpression::escape(baseName))); + QMap channelHeaders; // channelNum -> iprh path + + for (const QFileInfo &fi : dir.entryInfoList()) { + QRegularExpressionMatch match = reChannel.match(fi.fileName()); + if (match.hasMatch()) { + int chNum = match.captured(1).toInt(); + QString iprhPath = fi.absoluteFilePath(); + QString iprbPath = iprhPath; + iprbPath.replace(".iprh", ".iprb"); + if (QFile::exists(iprbPath)) { + channelHeaders.insert(chNum, iprhPath); + } else { + qDebug() << "Missing .iprb for" << iprhPath; + } + } + } + + if (channelHeaders.isEmpty()) { + qDebug() << "No multi-channel .iprh/.iprb found for baseName:" << baseName; + return false; + } + + const int channelCount = channelHeaders.size(); + qDebug() << "Impulse multi-channel: found" << channelCount << "channels for" << baseName; + + // 2. 解析第一个通道的头文件作为 master + auto it = channelHeaders.begin(); + GPRDataModel::Header masterHeader; + if (!parseIprhHeader(it.value(), masterHeader)) { + qDebug() << "Failed to parse master header:" << it.value(); + return false; + } + + QVector xOffsets; + QVector yOffsets; + xOffsets.reserve(channelCount); + yOffsets.reserve(channelCount); + + // 3. 验证各通道一致性并收集偏移量 + QVector channelTracesPerChannel; + channelTracesPerChannel.reserve(channelCount); + + double antennaSeparation = 0.0; + + for (auto cit = channelHeaders.begin(); cit != channelHeaders.end(); ++cit) { + int chNum = cit.key(); + QString iprhPath = cit.value(); + + GPRDataModel::Header chHeader; + if (!parseIprhHeader(iprhPath, chHeader)) { + qDebug() << "Failed to parse channel header:" << iprhPath; + return false; + } + + // 验证关键参数一致性 + if (chHeader.samplesPerTrace != masterHeader.samplesPerTrace) { + qDebug() << "Inconsistent SAMPLES across channels"; + return false; + } + if (qFuzzyCompare(chHeader.timeIntervalNs, masterHeader.timeIntervalNs) == false && + chHeader.timeIntervalNs > 0 && masterHeader.timeIntervalNs > 0) { + qDebug() << "Warning: Inconsistent TIME INTERVAL across channels"; + } + if (qFuzzyCompare(chHeader.distanceInc, masterHeader.distanceInc) == false && + chHeader.distanceInc > 0 && masterHeader.distanceInc > 0) { + qDebug() << "Warning: Inconsistent DISTANCE INTERVAL across channels"; + } + + // 收集单通道偏移量 + float xOff = 0.0f, yOff = 0.0f; + if (chHeader.rawParams.contains("CH_X_OFFSET")) { + xOff = chHeader.rawParams.value("CH_X_OFFSET").toFloat(); + } else if (chHeader.rawParams.contains("CH_OFFSET_X")) { + xOff = chHeader.rawParams.value("CH_OFFSET_X").toFloat(); + } else if (!chHeader.chXOffsets.isEmpty()) { + xOff = chHeader.chXOffsets.first(); + } + if (chHeader.rawParams.contains("CH_Y_OFFSET")) { + yOff = chHeader.rawParams.value("CH_Y_OFFSET").toFloat(); + } else if (chHeader.rawParams.contains("CH_OFFSET_Y")) { + yOff = chHeader.rawParams.value("CH_OFFSET_Y").toFloat(); + } else if (!chHeader.chYOffsets.isEmpty()) { + yOff = chHeader.chYOffsets.first(); + } + xOffsets.append(xOff); + yOffsets.append(yOff); + + if (chHeader.rawParams.contains("ANTENNA SEPARATION")) { + antennaSeparation = chHeader.rawParams.value("ANTENNA SEPARATION").toDouble(); + } + + // 从 .iprb 大小计算该通道的道数 + QString iprbPath = iprhPath; + iprbPath.replace(".iprh", ".iprb"); + QFile iprbFile(iprbPath); + int tracesInThisChannel = 0; + if (iprbFile.open(QIODevice::ReadOnly)) { + qint64 fileSize = iprbFile.size(); + qint64 traceBytes = static_cast(chHeader.samplesPerTrace) * sizeof(short); + if (traceBytes > 0) { + tracesInThisChannel = static_cast(fileSize / traceBytes); + } + iprbFile.close(); + } + channelTracesPerChannel.append(tracesInThisChannel); + } + + // 4. 以最小道数为准截断读取，容忍各通道微小差异 + int tracesPerChannel = channelTracesPerChannel.isEmpty() ? 0 : *std::min_element(channelTracesPerChannel.begin(), channelTracesPerChannel.end()); + for (int tc : channelTracesPerChannel) { + if (tc != tracesPerChannel) { + qDebug() << "Warning: Inconsistent trace count across channels. Using minimum:" << tracesPerChannel << "channel has:" << tc; + } + } + if (tracesPerChannel <= 0) { + qDebug() << "Error: No valid traces in any channel"; + return false; + } + + // 5. 如果偏移量全部为零，基于 ANTENNA SEPARATION 计算对称分布 + bool allZeroOffsets = true; + for (float v : xOffsets) { if (qFuzzyIsNull(v) == false) { allZeroOffsets = false; break; } } + if (allZeroOffsets && antennaSeparation > 1e-6 && channelCount > 1) { + double totalWidth = (channelCount - 1) * antennaSeparation; + double startX = -totalWidth / 2.0; + for (int i = 0; i < channelCount; ++i) { + xOffsets[i] = static_cast(startX + i * antennaSeparation); + } + qDebug() << "Computed symmetric X offsets from antenna separation:" << antennaSeparation; + } + + // 6. 组装合并后的 Header + model.header = masterHeader; + model.header.numberOfChannels = channelCount; + model.header.numTraces = tracesPerChannel * channelCount; + model.header.chXOffsets = xOffsets; + model.header.chYOffsets = yOffsets; + if (model.header.timeWindowNs <= 0.0 && model.header.timeIntervalNs > 0.0) { + model.header.timeWindowNs = model.header.samplesPerTrace * model.header.timeIntervalNs; + } + + // 7. 预分配 traces + model.traces.reserve(model.header.numTraces); + + // 8. 读取每个通道的 .iprb 到临时数组 + QVector> channelTraceArrays; + channelTraceArrays.resize(channelCount); + + int chIdx = 0; + for (auto cit = channelHeaders.begin(); cit != channelHeaders.end(); ++cit, ++chIdx) { + QString iprhPath = cit.value(); + QString iprbPath = iprhPath; + iprbPath.replace(".iprh", ".iprb"); + + QFile file(iprbPath); + if (!file.open(QIODevice::ReadOnly)) { + qDebug() << "Error: Cannot open" << iprbPath; + return false; + } + + const int samplesPerTrace = model.header.samplesPerTrace; + const qint64 singleTraceByteSize = samplesPerTrace * sizeof(short); + QByteArray traceBuffer; + traceBuffer.resize(singleTraceByteSize); + + channelTraceArrays[chIdx].reserve(tracesPerChannel); + + for (int t = 0; t < tracesPerChannel; ++t) { + if (file.atEnd()) break; + qint64 readBytes = file.read(traceBuffer.data(), traceBuffer.size()); + if (readBytes != traceBuffer.size()) break; + + RadarTrace oneTrace; + oneTrace.amplitudes.resize(samplesPerTrace); + const short* rawShortBuf = reinterpret_cast(traceBuffer.constData()); + for (int s = 0; s < samplesPerTrace; s++) { + oneTrace.amplitudes[s] = rawShortBuf[s]; + } + channelTraceArrays[chIdx].append(std::move(oneTrace)); + } + file.close(); + + if (channelTraceArrays[chIdx].size() != tracesPerChannel) { + qDebug() << "Warning: Channel" << cit.key() << "has fewer traces than expected"; + } + } + + // 9. 按 Mala Mira 格式交错合并 + model.channels = channelCount; + model.tracesPerChannel = tracesPerChannel; + for (int pos = 0; pos < tracesPerChannel; ++pos) { + for (int ch = 0; ch < channelCount; ++ch) { + if (pos < channelTraceArrays[ch].size()) { + RadarTrace trace = std::move(channelTraceArrays[ch][pos]); + trace.channelNumber = ch; + float xOff = (ch < xOffsets.size()) ? xOffsets[ch] : 0.0f; + float yOff = (ch < yOffsets.size()) ? yOffsets[ch] : 0.0f; + float lineDist = static_cast(model.header.startPosition + pos * model.header.distanceInc); + trace.position = QVector3D(xOff, lineDist - yOff, 0.0f); + model.traces.append(std::move(trace)); + } + } + } + + if (model.header.distanceInc > 1e-6) { + model.totalDistance = static_cast(model.tracesPerChannel * model.header.distanceInc); + } else { + model.totalDistance = static_cast(model.header.stopPosition - model.header.startPosition); + } + + qDebug() << "Impulse multi-channel load OK. Total traces:" << model.traces.size() + << "Channels:" << channelCount << "Traces/Channel:" << tracesPerChannel; + return !model.traces.isEmpty(); +} + +bool IprhParser::writeMalaFiles(const QString &radFilePath, + const QString &rd3FilePath, + const GPRDataModel &model, + const QString &sourceBaseName, + QString *errorMessage) +{ + QFile rd3File(rd3FilePath); + if (!rd3File.open(QIODevice::WriteOnly | QIODevice::Truncate)) { + if (errorMessage) { + *errorMessage = QStringLiteral("无法写入转换后的 RD3 文件：%1").arg(rd3FilePath); + } + return false; + } + + QDataStream out(&rd3File); + out.setByteOrder(QDataStream::LittleEndian); + const int samplesPerTrace = model.header.samplesPerTrace; + for (const RadarTrace &trace : model.traces) { + for (int s = 0; s < samplesPerTrace; ++s) { + const qint16 sample = static_cast(s < trace.amplitudes.size() ? trace.amplitudes[s] : 0); + out << sample; + } + } + rd3File.close(); + + QFile radFile(radFilePath); + if (!radFile.open(QIODevice::WriteOnly | QIODevice::Text | QIODevice::Truncate)) { + if (errorMessage) { + *errorMessage = QStringLiteral("无法写入转换后的 RAD 文件：%1").arg(radFilePath); + } + QFile::remove(rd3FilePath); + return false; + } + + QTextStream rad(&radFile); + rad << "# Converted from Impulse multi-channel survey: " << sourceBaseName << '\n'; + rad << "SAMPLES: " << model.header.samplesPerTrace << '\n'; + rad << "LAST TRACE: " << model.traces.size() << '\n'; + rad << "TIMEWINDOW: " << QLocale::c().toString(model.header.timeWindowNs, 'g', 12) << '\n'; + rad << "TIME INTERVAL: " << QLocale::c().toString(model.header.timeIntervalNs, 'g', 12) << '\n'; + rad << "DISTANCE INTERVAL: " << QLocale::c().toString(model.header.distanceInc, 'g', 12) << '\n'; + rad << "ANTENNAS: " << QLocale::c().toString(model.header.antennaFreq, 'g', 12) << '\n'; + if (!model.header.antennaType.isEmpty()) rad << "ANTENNA: " << model.header.antennaType << '\n'; + if (!model.header.date.isEmpty()) rad << "DATE: " << model.header.date << '\n'; + if (!model.header.timeStr.isEmpty()) rad << "TIME: " << model.header.timeStr << '\n'; + rad << "NUMBER_OF_CH: " << qMax(1, model.header.numberOfChannels) << '\n'; + rad << "CH_X_OFFSETS: " << formatFloatList(model.header.chXOffsets) << '\n'; + rad << "CH_Y_OFFSETS: " << formatFloatList(model.header.chYOffsets) << '\n'; + if (!model.header.units.isEmpty()) rad << "UNITS: " << model.header.units << '\n'; + rad << "START POSITION: " << QLocale::c().toString(model.header.startPosition, 'g', 12) << '\n'; + const double stopPosition = model.header.stopPosition > model.header.startPosition + ? model.header.stopPosition + : model.header.startPosition + model.tracesPerChannel * model.header.distanceInc; + rad << "STOP POSITION: " << QLocale::c().toString(stopPosition, 'g', 12) << '\n'; + radFile.close(); + + return true; +} + +QString IprhParser::formatFloatList(const QVector &values) +{ + QStringList parts; + parts.reserve(values.size()); + for (float value : values) { + parts.append(QLocale::c().toString(value, 'g', 10)); + } + return parts.join(' '); +} + +QString IprhParser::uniqueConvertedBasePath(const QString &dirPath, const QString &baseName) +{ + return QDir(dirPath).absoluteFilePath(baseName + QStringLiteral("_mala_converted")); +} + +double IprhParser::extractDouble(const QString &value) { + bool ok = false; + double res = value.toDouble(&ok); + return ok ? res : 0.0; +} + +int IprhParser::extractInt(const QString &value) { + bool ok = false; + int res = value.toInt(&ok); + return ok ? res : 0; +} diff --git a/external/gpr3dviewer/IprhParser.h b/external/gpr3dviewer/IprhParser.h new file mode 100644 index 0000000..7ae7add --- /dev/null +++ b/external/gpr3dviewer/IprhParser.h @@ -0,0 +1,41 @@ +#ifndef IPRHPARSER_H +#define IPRHPARSER_H + +#include "GPRDataModel.h" +#include +#include + +class IprhParser { +public: + // 主入口：传入 .iprh 文件路径，自动寻找同名的 .iprb 文件 + static bool loadFromIprh(const QString &iprhFilePath, GPRDataModel &model); + + // 多通道 Impulse 数据合并入口：一个文件夹下有多个 _A01.iprh/.iprb ... _A14.iprh/.iprb + static bool loadImpulseMultiChannel(const QString &dirPath, + const QString &baseName, + GPRDataModel &model); + + // 将多通道 Impulse 数据转换为本地 Mala Mira 风格 .rad/.rd3，再走现有 Mala Mira 解析流程 + static bool convertImpulseMultiChannelToMala(const QString &dirPath, + const QString &baseName, + QString *radFilePath, + QString *errorMessage = nullptr); + + static bool parseHeaderOnly(const QString &iprhFilePath, GPRDataModel::Header &header); + +private: + static bool parseIprhHeader(const QString &iprhFilePath, GPRDataModel::Header &header); + static bool loadIprbBinary(const QString &iprbFilePath, GPRDataModel &model); + static bool writeMalaFiles(const QString &radFilePath, + const QString &rd3FilePath, + const GPRDataModel &model, + const QString &sourceBaseName, + QString *errorMessage); + static QString formatFloatList(const QVector &values); + static QString uniqueConvertedBasePath(const QString &dirPath, const QString &baseName); + + static double extractDouble(const QString &value); + static int extractInt(const QString &value); +}; + +#endif // IPRHPARSER_H diff --git a/external/gpr3dviewer/PerformanceLogger.cpp b/external/gpr3dviewer/PerformanceLogger.cpp new file mode 100644 index 0000000..4c75e71 --- /dev/null +++ b/external/gpr3dviewer/PerformanceLogger.cpp @@ -0,0 +1,122 @@ +#include "PerformanceLogger.h" + +PerformanceLogger &PerformanceLogger::instance() +{ + static PerformanceLogger logger; + return logger; +} + +void PerformanceLogger::beginSession(const QString &sessionName) +{ + QMutexLocker locker(&m_mutex); + m_sessionName = sessionName; + m_sessionActive = true; +} + +void PerformanceLogger::endSession() +{ + QMutexLocker locker(&m_mutex); + m_sessionActive = false; +} + +void PerformanceLogger::log(const QString &category, const QString &name, qint64 elapsedMs) +{ + QMutexLocker locker(&m_mutex); + Record rec; + rec.category = category; + rec.name = name; + rec.elapsedMs = elapsedMs; + rec.timestamp = QDateTime::currentDateTime(); + m_records.append(rec); +} + +QVector PerformanceLogger::records() const +{ + QMutexLocker locker(&m_mutex); + return m_records; +} + +void PerformanceLogger::clear() +{ + QMutexLocker locker(&m_mutex); + m_records.clear(); +} + +QString PerformanceLogger::reportString() const +{ + QMutexLocker locker(&m_mutex); + QString report; + QTextStream ts(&report); + + ts << "==================================================\n"; + ts << "Performance Report"; + if (!m_sessionName.isEmpty()) + ts << " - " << m_sessionName; + ts << "\n"; + ts << "Generated: " << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss") << "\n"; + ts << "==================================================\n"; + + if (m_records.isEmpty()) { + ts << "No records.\n"; + return report; + } + + // Group by category + QString currentCategory; + qint64 categoryTotal = 0; + for (const auto &rec : m_records) { + if (rec.category != currentCategory) { + if (!currentCategory.isEmpty()) { + ts << " [Category Total: " << categoryTotal << " ms]\n\n"; + } + currentCategory = rec.category; + categoryTotal = 0; + ts << "[" << currentCategory << "]\n"; + } + ts << " " << rec.name << ": " << rec.elapsedMs << " ms\n"; + categoryTotal += rec.elapsedMs; + } + if (!currentCategory.isEmpty()) { + ts << " [Category Total: " << categoryTotal << " ms]\n"; + } + + ts << "\n==================================================\n"; + ts << "Grand Total: "; + qint64 grandTotal = 0; + for (const auto &rec : m_records) + grandTotal += rec.elapsedMs; + ts << grandTotal << " ms\n"; + ts << "==================================================\n"; + + return report; +} + +void PerformanceLogger::printReport() const +{ + const QString report = reportString(); + qDebug().noquote() << report; +} + +void PerformanceLogger::saveToFile(const QString &filePath) const +{ + QFile file(filePath); + if (!file.open(QIODevice::WriteOnly | QIODevice::Text | QIODevice::Append)) { + qDebug() << "PerformanceLogger: failed to open" << filePath; + return; + } + QTextStream ts(&file); + ts << reportString(); + file.close(); +} + +ScopedPerfTimer::ScopedPerfTimer(const QString &category, const QString &name) + : m_category(category), m_name(name) +{ + m_timer.start(); +} + +ScopedPerfTimer::~ScopedPerfTimer() +{ + const qint64 elapsed = m_timer.elapsed(); + PerformanceLogger::instance().log(m_category, m_name, elapsed); +} diff --git a/external/gpr3dviewer/PerformanceLogger.h b/external/gpr3dviewer/PerformanceLogger.h new file mode 100644 index 0000000..fbb237c --- /dev/null +++ b/external/gpr3dviewer/PerformanceLogger.h @@ -0,0 +1,60 @@ +#ifndef PERFORMANCELOGGER_H +#define PERFORMANCELOGGER_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +class PerformanceLogger +{ +public: + struct Record + { + QString category; + QString name; + qint64 elapsedMs = 0; + QDateTime timestamp; + }; + + static PerformanceLogger &instance(); + + void beginSession(const QString &sessionName = QString()); + void endSession(); + + void log(const QString &category, const QString &name, qint64 elapsedMs); + void printReport() const; + void saveToFile(const QString &filePath) const; + + QString reportString() const; + QVector records() const; + void clear(); + +private: + PerformanceLogger() = default; + mutable QMutex m_mutex; + QVector m_records; + QString m_sessionName; + bool m_sessionActive = false; +}; + +class ScopedPerfTimer +{ +public: + ScopedPerfTimer(const QString &category, const QString &name); + ~ScopedPerfTimer(); + +private: + QString m_category; + QString m_name; + QElapsedTimer m_timer; +}; + +#define SCOPED_PERF_TIMER(category, name) ScopedPerfTimer _perfTimer(category, name) + +#endif // PERFORMANCELOGGER_H diff --git a/external/gpr3dviewer/RadarProcessor.cpp b/external/gpr3dviewer/RadarProcessor.cpp new file mode 100644 index 0000000..c809429 --- /dev/null +++ b/external/gpr3dviewer/RadarProcessor.cpp @@ -0,0 +1,1148 @@ +/* + * RadarProcessor.cpp + * GPR探地雷达预处理算法实现源文件 + * 配套头文件 RadarProcessor.h + * 依赖库：kiss_fftr(快速傅里叶变换)、Qt容器、标准数学库 + * 整体流水线顺序：零时校正 → 零漂消除 → 背景压制 → TVG深度增益 → 带通滤波 → 剖面平滑 → 道内AGC → 道间均衡 → 希尔伯特包络/相位 + * 数据存储格式：short 16位整型振幅，全程运算double浮点，截断回short防止溢出 + */ +#include "RadarProcessor.h" +#include "GPRDataModel.h" +#include "PerformanceLogger.h" +#include +#include +#include +#include +#include +#include "kiss_fftr.h" +#include + +// 匿名命名空间：仅内部工具函数，对外不可见 +namespace { + +/** + * @brief 将浮点运算结果限制在short量程 [-32768,32767] + * @param value 运算后浮点振幅值 + * @return 钳位取整后的16位短整型 + */ +short clampToShort(double value) { + if (value > 32767.0) return 32767; + if (value < -32768.0) return -32768; + return static_cast(std::lround(value)); +} + +/** + * @brief 计算单道波形整体RMS均方根振幅 + * @param trace 单道雷达波形数据 + * @return 该道RMS能量值 + */ +double traceRms(const RadarTrace &trace) { + if (trace.amplitudes.isEmpty()) return 0.0; + double sumSquares = 0.0; + for (short amp : trace.amplitudes) { + double v = amp; + sumSquares += v * v; + } + return std::sqrt(sumSquares / trace.amplitudes.size()); +} + +/** + * @brief 全局整份数据所有采样点总RMS + * @param model 完整GPR多道数据集 + * @return 全局平均能量 + */ +double globalRms(const GPRDataModel &model) { + double sumSq = 0.0; + int cnt = 0; + for (const auto& tr : model.traces) { + for (short a : tr.amplitudes) { + double v = a; + sumSq += v * v; + cnt++; + } + } + if (cnt == 0) return 0.0; + return std::sqrt(sumSq / cnt); +} + +/** + * @brief 安全窗口尺寸限制，防止窗口超出数据长度 + * @param req 请求窗口大小 + * @param limit 最大允许长度 + * @return 合法区间[1,limit]内窗口值 + */ +int safeWindowSize(int req, int limit) { + if (limit <= 0) return 0; + return std::clamp(req, 1, limit); +} + +/** + * @brief 单道自动识别直达波零点起跳位置 + * @param trace 单道波形 + * @param frontWindow 仅在道前N个采样内搜索零点 + * @param sigmaMultiple 噪声标准差倍数起跳阈值 + * @return 零点对应的采样下标，-1代表识别失败 + * 逻辑：取道最前端小段计算基线噪声，使用正负极性对称的振幅+梯度条件判定起跳点 + */ +int detectSingleTraceZeroPoint(const RadarTrace &trace, int frontWindow, double sigmaMultiple) +{ + const int sampleCount = trace.amplitudes.size(); + const int searchEnd = std::min(frontWindow, sampleCount); + // 搜索区间太短无法判定基线噪声，直接返回失败 + if (searchEnd < 10) return -1; + + // 取搜索窗前1/4且不少于10点作为安静基线，避免固定10点对噪声估计过敏 + const int baselineLength = std::clamp(searchEnd / 4, 10, std::min(searchEnd, 40)); + double mean = 0.0; + for (int i = 0; i < baselineLength; ++i) { + mean += trace.amplitudes[i]; + } + mean /= baselineLength; + + // 计算基线方差、标准差 + double variance = 0.0; + for (int i = 0; i < baselineLength; ++i) { + const double diff = trace.amplitudes[i] - mean; + variance += diff * diff; + } + const double sigma = std::sqrt(variance / baselineLength); + const double threshold = std::max(1.0, sigmaMultiple * std::max(sigma, 1.0)); + const double gradientThreshold = std::max(1.0, 0.5 * threshold); + + // 逐点扫描梯度+振幅，正负极性直达波均可识别 + for (int sample = 1; sample < searchEnd; ++sample) { + const double current = trace.amplitudes[sample] - mean; + const double previous = trace.amplitudes[sample - 1] - mean; + const double gradient = current - previous; + if (std::abs(gradient) > gradientThreshold && std::abs(current) > threshold) { + return sample; + } + } + return -1; +} + +/** + * @brief 升余弦过渡带频域响应函数（带通滤波器窗函数） + * @param freq 当前频率点Hz + * @param low 通带下限Hz + * @param high 通带上限Hz + * @param transition 高低频过渡带宽Hz + * @return 频率增益系数[0,1] + * 作用：平缓滚降，消除矩形窗吉布斯振荡，滤波波形畸变更小 + */ +double bandpassResponse(double freq, double low, double high, double transition) +{ + // 升余弦渐入：低频侧从阻带到通带平滑抬升 + auto raisedCosineFadeIn = [](double t) { + return 0.5 * (1.0 - std::cos(M_PI * t)); + }; + // 升余弦渐出：高频侧从通带到阻带平滑衰减 + auto raisedCosineFadeOut = [](double t) { + return 0.5 * (1.0 + std::cos(M_PI * t)); + }; + + // 完全阻带，增益0 + if (freq < low - transition || freq > high + transition) return 0.0; + // 完全通带，增益1 + if (freq >= low + transition && freq <= high - transition) return 1.0; + + // 低频过渡区 + if (freq < low + transition) { + double t = (freq - (low - transition)) / (2.0 * transition); + t = std::clamp(t, 0.0, 1.0); + return raisedCosineFadeIn(t); + } + // 高频过渡区 + else { + double t = (freq - (high - transition)) / (2.0 * transition); + t = std::clamp(t, 0.0, 1.0); + return raisedCosineFadeOut(t); + } +} + +} // end 匿名命名空间 + +/** + * @brief 处理器构造函数，无资源初始化 + */ +RadarProcessor::RadarProcessor() +{ +} + +//===================================================================== +// 步骤0：零时自动/手动裁切 +//===================================================================== +/** + * @brief 切除每道零点前无效采样，统一道起始位置 + * @param input 原始输入数据 + * @param param 零时校正配置参数 + * @return 裁切后数据集 + * 自动模式：统计所有道零点中位数作为统一裁切长度；手动模式直接使用配置cutSamples + */ +GPRDataModel RadarProcessor::cutTimeZero(const GPRDataModel& input, const TimeZeroCutParams& param) +{ + SCOPED_PERF_TIMER("Processing", "cutTimeZero"); + // 开关未开启直接返回原数据 + if (!param.enable) return input; + GPRDataModel out = input; + int cutN = std::max(0, param.cutSamples); + + // 自动零点识别模式 + if (param.autoDetect) { + QVector zeroPoints; + zeroPoints.reserve(out.traces.size()); + const int nChannels = std::max(1, out.header.numberOfChannels); + QVector> channelZeroPoints(nChannels); + + // 逐道检测零点，失败值(-1)不参与统计，避免把全局中位数拉向0 + for (int ti = 0; ti < out.traces.size(); ++ti) { + const int zeroPoint = detectSingleTraceZeroPoint(out.traces[ti], param.frontSearchWindow, param.noiseSigmaMultiple); + if (zeroPoint <= 0) continue; + zeroPoints.append(zeroPoint); + channelZeroPoints[ti % nChannels].append(zeroPoint); + } + + QVector channelMedians; + channelMedians.reserve(nChannels); + for (QVector &points : channelZeroPoints) { + if (points.isEmpty()) continue; + std::sort(points.begin(), points.end()); + channelMedians.append(points[points.size() / 2]); + } + + if (!channelMedians.isEmpty()) { + // 多通道数据先按通道求中位数，再对通道中位数取中位数，避免通道间耦合差异互相污染 + std::sort(channelMedians.begin(), channelMedians.end()); + cutN = channelMedians[channelMedians.size() / 2]; + } else if (!zeroPoints.isEmpty()) { + std::sort(zeroPoints.begin(), zeroPoints.end()); + cutN = zeroPoints[zeroPoints.size() / 2]; + } + } + + // 无需裁切 + if (cutN <= 0) return out; + + // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach + for (auto &trace : out.traces) { + if (!trace.amplitudes.isEmpty()) + trace.amplitudes.data(); + } + + // 每道截断前面cutN个采样点 + RadarTrace *traces = out.traces.data(); + #pragma omp parallel for + for (int ti = 0; ti < out.traces.size(); ++ti) + { + RadarTrace &trace = traces[ti]; + int total = trace.amplitudes.size(); + if (cutN >= total) continue; + trace.amplitudes = trace.amplitudes.mid(cutN); + } + // 更新文件头每道采样总数 + if (!out.traces.isEmpty()) + { + out.header.samplesPerTrace = out.traces.first().amplitudes.size(); + } + return out; +} + +//===================================================================== +// 简易整道去直流（零漂DC模式底层实现） +//===================================================================== +/** + * @brief 扣除单道整体直流均值基线 + * @param input 输入数据 + * @param param 开关配置 + * @return 去直流后数据 + */ +GPRDataModel RadarProcessor::removeDCShift(const GPRDataModel &input, const DcShiftParams& param) +{ + SCOPED_PERF_TIMER("Processing", "removeDCShift"); + if (!param.enable) return input; + GPRDataModel output = input; + // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach + for (auto &trace : output.traces) { + if (!trace.amplitudes.isEmpty()) + trace.amplitudes.data(); + } + RadarTrace *traces = output.traces.data(); + #pragma omp parallel for + for (int i = 0; i < output.traces.size(); ++i) { + RadarTrace &trace = traces[i]; + if (trace.amplitudes.isEmpty()) continue; + // 计算单道整体均值 + double sum = 0; + for (auto amp : trace.amplitudes) sum += amp; + double mean = sum / trace.amplitudes.size(); + // 逐采样减去均值并钳位short + short *amps = trace.amplitudes.data(); + int sampleCount = trace.amplitudes.size(); + for (int s = 0; s < sampleCount; ++s) { + amps[s] = clampToShort(amps[s] - mean); + } + } + return output; +} + +//===================================================================== +// 零漂消除：DC全局均值 / Sliding滑动窗口基线 +//===================================================================== +/** + * @brief 消除波形基线漂移 + * @param input 输入数据 + * @param param 模式与窗口参数 + * @return 校正后数据 + * DC模式直接复用removeDCShift；滑动窗口用前缀和快速计算局部均值基线 + */ +GPRDataModel RadarProcessor::removeZeroDrift(const GPRDataModel &input, const ZeroDriftParams& param) +{ + SCOPED_PERF_TIMER("Processing", "removeZeroDrift"); + if (!param.enable) return input; + // 直流均值模式，调用简易去直流接口 + if (param.mode == ZeroDriftMode::DC) { + DcShiftParams dc; + dc.enable = true; + return removeDCShift(input, dc); + } + + GPRDataModel output = input; + // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach + for (auto &trace : output.traces) { + if (!trace.amplitudes.isEmpty()) + trace.amplitudes.data(); + } + RadarTrace *traces = output.traces.data(); + #pragma omp parallel for + for (int ti = 0; ti < output.traces.size(); ++ti) { + RadarTrace &trace = traces[ti]; + const QVector source = trace.amplitudes; + const int sampleCount = source.size(); + if (sampleCount == 0) continue; + + const int window = safeWindowSize(param.slidingWindowSamples, sampleCount); + const int halfWindow = window / 2; + // 前缀和数组，O(1)区间求和加速滑动窗口均值 + QVector prefix(sampleCount + 1, 0.0); + for (int sample = 0; sample < sampleCount; ++sample) { + prefix[sample + 1] = prefix[sample] + source[sample]; + } + + // 逐采样计算窗口内局部基线并扣除 + short *amps = trace.amplitudes.data(); + for (int sample = 0; sample < sampleCount; ++sample) { + const int start = std::max(0, sample - halfWindow); + const int end = std::min(sampleCount, start + window); + const int adjStart = std::max(0, end - window); + const int count = end - adjStart; + const double localMean = count > 0 ? (prefix[end] - prefix[adjStart]) / count : 0.0; + amps[sample] = clampToShort(source[sample] - localMean); + } + } + return output; +} + +//===================================================================== +// 全局统一倍率增益 +//===================================================================== +/** + * @brief 整份数据统一放大缩小倍数 + * @param input 输入数据 + * @param param 增益系数与开关 + * @return 缩放后数据 + */ +GPRDataModel RadarProcessor::applyGain(const GPRDataModel &input, const GlobalGainParams& param) +{ + SCOPED_PERF_TIMER("Processing", "applyGain"); + if (!param.enable) return input; + GPRDataModel output = input; + double f = param.factor; + // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach + for (auto &trace : output.traces) { + if (!trace.amplitudes.isEmpty()) + trace.amplitudes.data(); + } + RadarTrace *traces = output.traces.data(); + #pragma omp parallel for + for (int ti = 0; ti < output.traces.size(); ++ti) { + RadarTrace &tr = traces[ti]; + short *amps = tr.amplitudes.data(); + int sampleCount = tr.amplitudes.size(); + for (int s = 0; s < sampleCount; ++s) { + amps[s] = clampToShort(amps[s] * f); + } + } + return output; +} + +//===================================================================== +// TVG 球面扩散+介质吸收深度增益补偿 +//===================================================================== +/** + * @brief 随深度自动补偿电磁波衰减 + * @param input 原始数据 + * @param params 速度、吸收系数、最大增益限制等配置 + * @return 增益补偿后剖面 + * 公式：增益 = 球面扩散项 × 指数吸收项，上限maxGain防止噪声爆炸放大 + */ +GPRDataModel RadarProcessor::sphericalTvg(const GPRDataModel &input, const SphericalTvgParams ¶ms) +{ + SCOPED_PERF_TIMER("Processing", "sphericalTvg"); + if (!params.enable) return input; + GPRDataModel output = input; + if (output.traces.isEmpty() || output.header.samplesPerTrace <= 0 || output.header.timeWindowNs <= 0.0) return output; + + // 波速优先级：参数指定 > 文件头读取 > 兜底0.1 m/ns + const double velocity = params.velocityMPerNs > 0.0 ? params.velocityMPerNs : (output.header.waveVelocity > 0.0 ? output.header.waveVelocity : 0.1); + const double referenceDepth = std::max(params.referenceDepthM, std::numeric_limits::epsilon()); + const double exponent = std::max(params.exponent, 0.0); + const double maxGain = params.maxGain > 0.0 ? params.maxGain : 1.0; + // 单采样时间间隔ns + const double dtNs = output.header.timeWindowNs / output.header.samplesPerTrace; + + // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach + for (auto &trace : output.traces) { + if (!trace.amplitudes.isEmpty()) + trace.amplitudes.data(); + } + RadarTrace *traces = output.traces.data(); + #pragma omp parallel for + for (int ti = 0; ti < output.traces.size(); ++ti) { + RadarTrace &trace = traces[ti]; + short *amps = trace.amplitudes.data(); + const int sampleCount = trace.amplitudes.size(); + for (int sample = 0; sample < sampleCount; ++sample) { + const double timeNs = sample * dtNs; + // 单程深度 = 往返时间 × 波速 / 2 + const double depthM = timeNs * velocity / 2.0; + double gain = 1.0; + // 球面扩散补偿 (r/r0)^n + if (params.enableSpherical) { + gain *= std::pow(std::max(depthM, referenceDepth) / referenceDepth, exponent); + } + // 介质指数吸收补偿 exp(β·t) + if (params.enableAbsorption && params.absorptionBeta > 0.0) { + gain *= std::exp(params.absorptionBeta * timeNs); + } + // 增益封顶 + gain = std::min(gain, maxGain); + amps[sample] = clampToShort(amps[sample] * gain); + } + } + return output; +} + +//===================================================================== +// 二维剖面平滑：垂直(深度) + 水平(测线道)均值平滑 +//===================================================================== +/** + * @brief 时空二维滑动均值降噪 + * @param input 输入数据集 + * @param params 窗口半宽、横竖平滑开关 + * @return 平滑后剖面 + * 先深度方向逐道平滑，再道方向同通道横向平滑；横向使用原始未平滑数据做基准避免二次模糊 + */ +GPRDataModel RadarProcessor::profileSmooth(const GPRDataModel &input, const ProfileSmoothParams ¶ms) +{ + SCOPED_PERF_TIMER("Processing", "profileSmooth"); + if (!params.enable) return input; + GPRDataModel output = input; + if (output.traces.isEmpty()) return output; + + const int nChannels = output.header.numberOfChannels > 0 ? output.header.numberOfChannels : 1; + const int tracesPerChannel = output.getTracesPerChannel(); + const int winHalf = std::max(0, params.smoothWindow); + if (winHalf == 0) return output; + + // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach + for (auto &trace : output.traces) { + if (!trace.amplitudes.isEmpty()) + trace.amplitudes.data(); + } + + // 第一步：垂直深度方向单道内平滑 + if (params.verticalSmooth) { + RadarTrace *traces = output.traces.data(); + #pragma omp parallel for + for (int ti = 0; ti < output.traces.size(); ++ti) { + RadarTrace &trace = traces[ti]; + const QVector source = trace.amplitudes; + const int sampleCount = source.size(); + short *amps = trace.amplitudes.data(); + for (int sample = 0; sample < sampleCount; ++sample) { + // 窗口左右边界钳位 + const int start = std::max(0, sample - winHalf); + const int end = std::min(sampleCount - 1, sample + winHalf); + double sum = 0.0; + for (int idx = start; idx <= end; ++idx) { + sum += source[idx]; + } + amps[sample] = clampToShort(sum / (end - start + 1)); + } + } + } + + // 第二步：水平测线方向同通道多道平滑，使用原始模型做参考源 + if (params.horizontalSmooth && tracesPerChannel > 1) { + const GPRDataModel sourceModel = output; + RadarTrace *traces = output.traces.data(); + #pragma omp parallel for + for (int globalIndex = 0; globalIndex < output.traces.size(); ++globalIndex) { + int channel = globalIndex % nChannels; + int traceInChannel = globalIndex / nChannels; + + RadarTrace &targetTrace = traces[globalIndex]; + const int sampleCount = targetTrace.amplitudes.size(); + const int startTrace = std::max(0, traceInChannel - winHalf); + const int endTrace = std::min(tracesPerChannel - 1, traceInChannel + winHalf); + // 每个采样点横向多道平均 + for (int sample = 0; sample < sampleCount; ++sample) { + double sum = 0.0; + int count = 0; + for (int neighbor = startTrace; neighbor <= endTrace; ++neighbor) { + const int neighborIndex = sourceModel.getTraceIndex(channel, neighbor); + if (neighborIndex < 0 || neighborIndex >= sourceModel.traces.size()) continue; + const QVector &s = sourceModel.traces[neighborIndex].amplitudes; + if (sample >= amps.size()) continue; + sum += amps[sample]; + ++count; + } + if (count > 0) { + short *tgtAmps = targetTrace.amplitudes.data(); + tgtAmps[sample] = clampToShort(sum / count); + } + } + } + } + return output; +} + +//===================================================================== +// 道内滑动窗口AGC振幅均衡 +//===================================================================== +/** + * @brief 单道内部深浅振幅自适应均衡 + * @param input 原始数据 + * @param params 窗口大小、目标RMS、最大增益限制 + * @return AGC均衡后数据 + * 前缀平方和快速计算局部RMS，增益=目标RMS/局部RMS，封顶maxGain防止噪声放大 + */ +GPRDataModel RadarProcessor::traceInnerAgc(const GPRDataModel &input, const TraceInnerAgcParams ¶ms) +{ + SCOPED_PERF_TIMER("Processing", "traceInnerAgc"); + if (!params.enable) return input; + GPRDataModel output = input; + if (output.traces.isEmpty()) return output; + + // 目标RMS：0则自动取整份数据全局RMS + const double targetRms = params.targetRms > 0.0 ? params.targetRms : globalRms(input); + if (targetRms <= 0.0) return output; + const double epsilon = params.epsilon > 0.0 ? params.epsilon : 1.0; + const double maxGain = params.maxGain > 0.0 ? params.maxGain : 1.0; + + // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach + for (auto &trace : output.traces) { + if (!trace.amplitudes.isEmpty()) + trace.amplitudes.data(); + } + const RadarTrace *inputTraces = input.traces.data(); + RadarTrace *outputTraces = output.traces.data(); + #pragma omp parallel for + for (int traceIndex = 0; traceIndex < input.traces.size(); ++traceIndex) { + const RadarTrace &sourceTrace = inputTraces[traceIndex]; + RadarTrace &targetTrace = outputTraces[traceIndex]; + const int sampleCount = sourceTrace.amplitudes.size(); + if (sampleCount == 0) continue; + + const int window = safeWindowSize(params.windowSamples, sampleCount); + const int halfWindow = window / 2; + // 平方前缀和数组，快速区间RMS计算 + QVector prefixSquares(sampleCount + 1, 0.0); + const short *srcAmps = sourceTrace.amplitudes.data(); + for (int s = 0; s < sampleCount; s++) { + double v = srcAmps[s]; + prefixSquares[s+1] = prefixSquares[s] + v * v; + } + // 逐采样计算局部增益并缩放振幅 + short *tgtAmps = targetTrace.amplitudes.data(); + for (int s = 0; s < sampleCount; s++) { + int start = std::max(0, s - halfWindow); + int end = std::min(sampleCount, start + window); + int adjStart = std::max(0, end - window); + int cnt = end - adjStart; + double sumSq = prefixSquares[end] - prefixSquares[adjStart]; + double localRms = cnt > 0 ? std::sqrt(sumSq / cnt) : 0.0; + // 防除零极小值epsilon + double gain = std::min(targetRms / std::max(localRms, epsilon), maxGain); + tgtAmps[s] = clampToShort(srcAmps[s] * gain); + } + } + return output; +} + +//===================================================================== +// 道间AGC：整条测线道与道振幅均衡 +//===================================================================== +/** + * @brief 横向测线多道之间整体能量拉平 + * @param input 输入数据 + * @param params 全局/局部滑动窗口模式 + * @return 均衡后剖面 + * Global：全部道共用一个基准RMS；Local：相邻窗口道平均RMS做基准 + */ +GPRDataModel RadarProcessor::traceToTraceEqualization(const GPRDataModel &input, const TraceToTraceAgcParams ¶ms) +{ + SCOPED_PERF_TIMER("Processing", "traceToTraceEqualization"); + if (!params.enable) return input; + GPRDataModel output = input; + int traceCnt = input.traces.size(); + if (traceCnt == 0) return output; + + const int nChannels = std::max(1, input.header.numberOfChannels); + const int tracesPerChannel = input.getTracesPerChannel(); + + // 预计算每道整体RMS能量 + QVector trRms(traceCnt, 0.0); + for (int i = 0; i < traceCnt; i++) trRms[i] = traceRms(input.traces[i]); + double targetRms = params.targetRms > 0 ? params.targetRms : globalRms(input); + if (targetRms <= 0) return output; + + double eps = params.epsilon > 0 ? params.epsilon : 1.0; + double maxG = params.maxGain > 0 ? params.maxGain : 1.0; + int win = safeWindowSize(params.horizontalWindowTraces, std::max(1, tracesPerChannel)); + int halfWin = win / 2; + + QVector> channelPrefixSq(nChannels); + if (params.mode == TraceToTraceMode::Local && tracesPerChannel > 0) { + for (int ch = 0; ch < nChannels; ++ch) { + channelPrefixSq[ch].resize(tracesPerChannel + 1); + channelPrefixSq[ch][0] = 0.0; + for (int tr = 0; tr < tracesPerChannel; ++tr) { + const int idx = input.getTraceIndex(ch, tr); + const double rms = (idx >= 0 && idx < traceCnt) ? trRms[idx] : 0.0; + channelPrefixSq[ch][tr + 1] = channelPrefixSq[ch][tr] + rms * rms; + } + } + } + + // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach + for (auto &trace : output.traces) { + if (!trace.amplitudes.isEmpty()) + trace.amplitudes.data(); + } + RadarTrace *traces = output.traces.data(); + #pragma omp parallel for + for (int ti = 0; ti < traceCnt; ti++) + { + double base = trRms[ti]; + // 局部模式只沿同一通道的测线方向取滑动窗口RMS，避免多通道交错数据互相污染 + if (params.mode == TraceToTraceMode::Local && tracesPerChannel > 0) + { + const int channel = ti % nChannels; + const int traceInChannel = ti / nChannels; + if (channel >= 0 && channel < channelPrefixSq.size() && traceInChannel >= 0 && traceInChannel < tracesPerChannel) { + int st = std::max(0, traceInChannel - halfWin); + int ed = std::min(tracesPerChannel, st + win); + int adjSt = std::max(0, ed - win); + int c = ed - adjSt; + double sq = channelPrefixSq[channel][ed] - channelPrefixSq[channel][adjSt]; + base = c > 0 ? std::sqrt(sq / c) : base; + } + } + // 计算整道统一增益 + double g = std::min(targetRms / std::max(base, eps), maxG); + // 整道所有采样同倍率缩放 + short *amps = traces[ti].amplitudes.data(); + int sampleCount = traces[ti].amplitudes.size(); + for (int s = 0; s < sampleCount; ++s) + { + amps[s] = clampToShort(amps[s] * g); + } + } + return output; +} + +//===================================================================== +// 希尔伯特变换：振幅包络 / 正交虚部 +//===================================================================== +/** + * @brief 时域直接希尔伯特变换（简易实现，适合中小采样长度） + * @param input 输入波形 + * @param params 输出包络/正交分量开关 + * @return 变换后数据 + * 公式：虚部h[n] = Σ(2/(π·(n-k)))·x[k]；包络=√(实部²+虚部²) + */ +GPRDataModel RadarProcessor::hilbertTransform(const GPRDataModel &input, const HilbertParams ¶ms) +{ + SCOPED_PERF_TIMER("Processing", "hilbertTransform"); + if (!params.enable) return input; + GPRDataModel output = input; + if (output.traces.isEmpty()) return output; + + // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach + for (auto &trace : output.traces) { + if (!trace.amplitudes.isEmpty()) + trace.amplitudes.data(); + } + RadarTrace *traces = output.traces.data(); + #pragma omp parallel for + for (int ti = 0; ti < output.traces.size(); ++ti) { + RadarTrace &trace = traces[ti]; + const QVector source = trace.amplitudes; + const int sampleCount = source.size(); + if (sampleCount < 3) continue; + + QVector hilbert(sampleCount, 0.0); + // 时域卷积计算希尔伯特虚部 + for (int n = 0; n < sampleCount; ++n) { + double sum = 0.0; + for (int k = 0; k < sampleCount; ++k) { + const int diff = n - k; + // 差分0、偶数项系数为0，只累加奇数差分 + if (diff == 0 || diff % 2 == 0) continue; + sum += (2.0 / (M_PI * diff)) * source[k]; + } + hilbert[n] = sum; + } + + // 输出模式分支 + short *amps = trace.amplitudes.data(); + for (int sample = 0; sample < sampleCount; ++sample) { + if (params.computeEnvelope) { + // 成像常用：振幅包络 + const double real = source[sample]; + amps[sample] = clampToShort(std::sqrt(real * real + hilbert[sample] * hilbert[sample])); + } else { + // 相位分析：仅输出正交虚部 + amps[sample] = clampToShort(hilbert[sample]); + } + } + } + return output; +} + +//===================================================================== +// FFT频域带通滤波（kiss_fftr快速实数FFT） +//===================================================================== +/** + * @brief 频域升余弦窗带通滤波 + * @param input 原始数据 + * @param params 自动/手动频带配置 + * @return 滤波后波形 + * 流程：时域→FFT频域乘增益响应→逆FFT回时域；补零2倍长度避免循环卷积混叠 + */ +GPRDataModel RadarProcessor::bandpassFilter(const GPRDataModel &input, const BandpassParams ¶ms) +{ + SCOPED_PERF_TIMER("Processing", "bandpassFilter"); + if (!params.enable) return input; + GPRDataModel out = input; + if (out.traces.isEmpty() || out.header.samplesPerTrace <= 0) return out; + + const double dtNs = out.header.timeWindowNs / out.header.samplesPerTrace; + if (dtNs <= 0.0) return out; + const double fs = 1e9 / dtNs; + + double lowFreq = params.lowFreqHz; + double highFreq = params.highFreqHz; + if (params.autoFreq && params.antennaFreqMHz > 0.0) { + const double fcenter = params.antennaFreqMHz * 1e6; + lowFreq = fcenter * 0.1; + highFreq = fcenter * 2.5; + } + lowFreq = std::max(lowFreq, 0.0); + highFreq = std::min(highFreq, fs * 0.5); + if (lowFreq >= highFreq) return out; + + const double transition = std::max(50.0, (highFreq - lowFreq) * 0.1); + const int sampleCount = out.header.samplesPerTrace; + + int nfft = 1; + while (nfft < sampleCount * 2) nfft <<= 1; + + std::vector filterResponse(nfft / 2 + 1); + for (int k = 0; k <= nfft / 2; ++k) { + const double freq = k * fs / nfft; + filterResponse[k] = static_cast(bandpassResponse(freq, lowFreq, highFreq, transition)); + } + + for (auto &trace : out.traces) { + if (!trace.amplitudes.isEmpty()) + trace.amplitudes.data(); + } + + #pragma omp parallel + { + kiss_fftr_cfg fwd = kiss_fftr_alloc(nfft, 0, nullptr, nullptr); + kiss_fftr_cfg inv = kiss_fftr_alloc(nfft, 1, nullptr, nullptr); + std::vector timedata(nfft, 0.0f); + std::vector freqdata(nfft / 2 + 1); + + RadarTrace *traces = out.traces.data(); + if (fwd && inv) { + #pragma omp for + for (int ti = 0; ti < out.traces.size(); ++ti) { + RadarTrace &trace = traces[ti]; + if (trace.amplitudes.size() < 4) continue; + + short *amps = trace.amplitudes.data(); + for (int i = 0; i < sampleCount; ++i) { + timedata[i] = static_cast(amps[i]); + } + for (int i = sampleCount; i < nfft; ++i) { + timedata[i] = 0.0f; + } + + kiss_fftr(fwd, timedata.data(), freqdata.data()); + + for (int k = 0; k <= nfft / 2; ++k) { + freqdata[k].r *= filterResponse[k]; + freqdata[k].i *= filterResponse[k]; + } + + kiss_fftri(inv, freqdata.data(), timedata.data()); + + const double scale = 1.0 / nfft; + for (int i = 0; i < sampleCount; ++i) { + amps[i] = clampToShort(timedata[i] * scale); + } + } + } + + if (fwd) kiss_fftr_free(fwd); + if (inv) kiss_fftr_free(inv); + } + return out; +} + +//===================================================================== +// 背景杂波去除：均值法 / 奇异值过滤法 +//===================================================================== +/** + * @brief 多道平均背景相减，压制固定耦合、地层静态杂波 + * @param input 原始剖面 + * @param params 模式、平均窗口道数、奇异阈值 + * @return 去除背景后数据 + * MeanAverage：全部道参与背景平均；SingularityFilter：过滤高能异常道再平均，保护空洞/管线强反射 + */ +GPRDataModel RadarProcessor::backgroundRemoval(const GPRDataModel &input, const BackgroundParams ¶ms) +{ + SCOPED_PERF_TIMER("Processing", "backgroundRemoval"); + if (!params.enable) return input; + GPRDataModel out = input; + const int traceTotal = out.traces.size(); + if (traceTotal < 3) return out; + + const int nChannels = std::max(1, out.header.numberOfChannels); + const int tracesPerChannel = out.getTracesPerChannel(); + if (tracesPerChannel < 3) return out; + + QVector> validTraceMask(nChannels); + for (int ch = 0; ch < nChannels; ++ch) { + validTraceMask[ch].fill(true, tracesPerChannel); + } + + // 奇异滤波模式：每个通道内部筛掉能量过高的异常道，只用平稳道算背景 + if (params.mode == BackgroundMode::SingularityFilter) { + for (int ch = 0; ch < nChannels; ++ch) { + QVector energies; + energies.reserve(tracesPerChannel); + for (int tr = 0; tr < tracesPerChannel; ++tr) { + const int idx = input.getTraceIndex(ch, tr); + if (idx >= 0 && idx < traceTotal) { + energies.append(traceRms(input.traces[idx])); + } + } + if (energies.isEmpty()) continue; + + QVector sortedEnergies = energies; + std::sort(sortedEnergies.begin(), sortedEnergies.end()); + const double median = sortedEnergies[sortedEnergies.size() / 2]; + const double limit = median * std::max(params.singularityThreshold, 1.0); + + for (int tr = 0; tr < tracesPerChannel; ++tr) { + const int idx = input.getTraceIndex(ch, tr); + validTraceMask[ch][tr] = idx >= 0 && idx < traceTotal && traceRms(input.traces[idx]) <= limit; + } + } + } + + // 预先detach所有trace的amplitudes，避免OpenMP多线程竞争隐式共享容器的detach + for (auto &trace : out.traces) { + if (!trace.amplitudes.isEmpty()) + trace.amplitudes.data(); + } + + // 滑动半道窗口只沿同一通道的测线方向取背景，避免多通道交错数据互相污染 + const int halfWindow = std::max(1, params.averageTraceCount) / 2; + RadarTrace *outTraces = out.traces.data(); + const RadarTrace *inTraces = input.traces.data(); + #pragma omp parallel for + for (int ti = 0; ti < traceTotal; ++ti) + { + const int channel = ti % nChannels; + const int traceInChannel = ti / nChannels; + if (channel < 0 || channel >= validTraceMask.size() || traceInChannel < 0 || traceInChannel >= tracesPerChannel) continue; + + const int left = std::max(0, traceInChannel - halfWindow); + const int right = std::min(tracesPerChannel - 1, traceInChannel + halfWindow); + const int sampleCount = outTraces[ti].amplitudes.size(); + QVector background(sampleCount, 0.0); + int backgroundCount = 0; + + for (int tr = left; tr <= right; ++tr) { + if (!validTraceMask[channel][tr]) continue; + const int idx = input.getTraceIndex(channel, tr); + if (idx < 0 || idx >= traceTotal) continue; + const RadarTrace &bgTrace = inTraces[idx]; + const short *bgAmps = bgTrace.amplitudes.data(); + int bgSize = bgTrace.amplitudes.size(); + for (int sample = 0; sample < sampleCount && sample < bgSize; ++sample) { + background[sample] += bgAmps[sample]; + } + ++backgroundCount; + } + if (backgroundCount == 0) continue; + + RadarTrace ¤tTrace = outTraces[ti]; + const QVector &source = inTraces[ti].amplitudes; + short *curAmps = currentTrace.amplitudes.data(); + const short *srcAmps = source.data(); + int srcSize = source.size(); + for (int sample = 0; sample < sampleCount && sample < srcSize; ++sample) + { + const double averageBackground = background[sample] / backgroundCount; + curAmps[sample] = clampToShort(srcAmps[sample] - averageBackground); + } + } + return out; +} + +//===================================================================== +// Kirchhoff偏移：x-t域双曲线求和 +//===================================================================== +/** + * @brief Kirchhoff偏移在x-t域实现 + * @param input 输入剖面数据 + * @param params 求和宽度（两侧道数）与波速参数 + * @return 偏移后数据 + * 算法流程： + * 1. 对原始剖面逐道时间方向微分 + * 2. 对每个输出点(t0, x0)沿双曲线 t = sqrt(t0^2 + (dx*n/v)^2) 求和 + * 3. 加权因子采用1/t几何扩散补偿 + * 4. sumWidth < 128时预计算双曲线查表，避免重复sqrt运算 + */ +GPRDataModel RadarProcessor::kirchhoffMigration(const GPRDataModel &input, const MigrationParams ¶ms) +{ + SCOPED_PERF_TIMER("Processing", "kirchhoffMigration"); + if (!params.enable) return input; + GPRDataModel out = input; + int traceTotal = out.traces.size(); + if (traceTotal < 3) return out; + int sampleCount = out.header.samplesPerTrace; + if (sampleCount < 3) return out; + + const double dtNs = out.header.timeWindowNs / sampleCount; + double dx = out.header.distanceInc; + if (dx <= 0.0) dx = 0.01; + double v = params.velocityMPerNs > 0.0 ? params.velocityMPerNs : out.header.waveVelocity; + if (v <= 0.0) v = 0.1; + + const int sumWidth = std::max(1, params.sumWidth); + + // 1. 逐道时间方向微分，结果存为double浮点数组 + QVector> diffData(traceTotal); + for (int ti = 0; ti < traceTotal; ++ti) { + diffData[ti].resize(sampleCount); + const QVector &src = out.traces[ti].amplitudes; + if (sampleCount > 0) diffData[ti][0] = 0.0; + for (int s = 1; s < sampleCount; ++s) { + diffData[ti][s] = static_cast(src[s]) - static_cast(src[s - 1]); + } + } + + // 2. 当求和宽度小于128时，预计算双曲线采样位置查表（只进行一次） + QVector> hyperbolaTable; + if (sumWidth < 128) { + hyperbolaTable.resize(sampleCount); + for (int t0 = 0; t0 < sampleCount; ++t0) { + hyperbolaTable[t0].resize(sumWidth * 2 + 1); + const double t0Ns = t0 * dtNs; + for (int k = -sumWidth; k <= sumWidth; ++k) { + const double x = k * dx; + const double tNs = std::sqrt(t0Ns * t0Ns + (x / v) * (x / v)); + hyperbolaTable[t0][k + sumWidth] = tNs / dtNs; + } + } + } + + // 预先detach所有trace的amplitudes + for (auto &trace : out.traces) { + if (!trace.amplitudes.isEmpty()) + trace.amplitudes.data(); + } + + // 3. 对每个输出点沿双曲线加权求和 + RadarTrace *traces = out.traces.data(); + #pragma omp parallel for + for (int ti = 0; ti < traceTotal; ++ti) { + QVector result(sampleCount, 0); + for (int s = 0; s < sampleCount; ++s) { + double sum = 0.0; + int validCount = 0; + for (int k = -sumWidth; k <= sumWidth; ++k) { + const int srcTrace = ti + k; + if (srcTrace < 0 || srcTrace >= traceTotal) continue; + + double srcSampleD; + if (sumWidth < 128 && !hyperbolaTable.isEmpty()) { + srcSampleD = hyperbolaTable[s][k + sumWidth]; + } else { + const double t0Ns = s * dtNs; + const double x = k * dx; + const double tNs = std::sqrt(t0Ns * t0Ns + (x / v) * (x / v)); + srcSampleD = tNs / dtNs; + } + + int srcSample = static_cast(srcSampleD); + double frac = srcSampleD - srcSample; + if (srcSample < 0) { + srcSample = 0; + frac = 0.0; + } + if (srcSample >= sampleCount - 1) { + srcSample = sampleCount - 1; + frac = 0.0; + } + + double amp = diffData[srcTrace][srcSample]; + if (frac > 0.0 && srcSample < sampleCount - 1) { + amp += frac * (diffData[srcTrace][srcSample + 1] - diffData[srcTrace][srcSample]); + } + + // 1/t 几何扩散加权（t>0时） + const double tNs = s * dtNs; + if (tNs > 1e-12) { + amp /= tNs; + } + + sum += amp; + ++validCount; + } + if (validCount > 0) { + result[s] = clampToShort(sum); + } + } + traces[ti].amplitudes = std::move(result); + } + + return out; +} + +//===================================================================== +// 流水线总调度入口：按步骤队列顺序串行执行所有开启算法 +//===================================================================== +/** + * @brief 批量串行执行整套预处理流水线 + * @param rawData 未处理原始GPR数据 + * @param pipeline 有序步骤配置容器 + * @return 逐步骤处理完成后的成品数据 + */ +GPRDataModel RadarProcessor::runPipeline(const GPRDataModel& rawData, const ProcPipeline& pipeline) +{ + return runPipeline(rawData, pipeline, nullptr); +} + +GPRDataModel RadarProcessor::runPipeline(const GPRDataModel& rawData, + const ProcPipeline& pipeline, + const PipelineRuntimeContext *context) +{ + SCOPED_PERF_TIMER("Processing", "runPipeline(Total)"); + auto canceled = [context]() { + return context && context->isCanceled && context->isCanceled(); + }; + auto stepName = [](ProcStepType type) -> QString { + switch (type) { + case ProcStepType::StepTimeZeroCut: return QStringLiteral("Set Time Zero"); + case ProcStepType::StepZeroDrift: return QStringLiteral("Zero Drift"); + case ProcStepType::StepRemoveDC: return QStringLiteral("Remove DC"); + case ProcStepType::StepBackgroundRemove: return QStringLiteral("Background Removal"); + case ProcStepType::StepSphericalTVG: return QStringLiteral("Spherical TVG"); + case ProcStepType::StepBandpassFilter: return QStringLiteral("Bandpass Filter"); + case ProcStepType::StepProfileSmooth: return QStringLiteral("Profile Smooth"); + case ProcStepType::StepInnerAGC: return QStringLiteral("Trace Inner AGC"); + case ProcStepType::StepTraceToTraceAGC: return QStringLiteral("Trace-to-Trace AGC"); + case ProcStepType::StepHilbertTransform: return QStringLiteral("Hilbert Transform"); + case ProcStepType::StepGlobalGain: return QStringLiteral("Global Gain"); + case ProcStepType::StepMigration: return QStringLiteral("Migration"); + default: return QStringLiteral("Unknown Step"); + } + }; + + GPRDataModel current = rawData; + // 处理流水线只操作 traces，volumeData 不参与任何算法步骤 + // 提前清除可显著减少各步骤内部 GPRDataModel out = input 时的深拷贝开销 + current.volumeData.clear(); + // 严格按照steps数组存储顺序依次执行 + const int stepCount = pipeline.steps.size(); + for (int i = 0; i < pipeline.steps.size(); ++i) + { + if (canceled()) + break; + + const auto& unit = pipeline.steps[i]; + qDebug() << "[RadarProcessor] Step" << (i + 1) << "/" << stepCount << "starting:" << stepName(unit.type) + << "| traces=" << current.traces.size() << "samples=" << current.header.samplesPerTrace; + + switch (unit.type) + { + case ProcStepType::StepTimeZeroCut: + current = cutTimeZero(current, unit.zeroCut); + break; + case ProcStepType::StepZeroDrift: + current = removeZeroDrift(current, unit.zeroDrift); + break; + case ProcStepType::StepSphericalTVG: + current = sphericalTvg(current, unit.tvg); + break; + case ProcStepType::StepBandpassFilter: + current = bandpassFilter(current, unit.band); + break; + case ProcStepType::StepProfileSmooth: + current = profileSmooth(current, unit.smooth); + break; + case ProcStepType::StepBackgroundRemove: + current = backgroundRemoval(current, unit.bg); + break; + case ProcStepType::StepRemoveDC: + current = removeDCShift(current, unit.dc); + break; + case ProcStepType::StepInnerAGC: + current = traceInnerAgc(current, unit.innerAgc); + break; + case ProcStepType::StepTraceToTraceAGC: + current = traceToTraceEqualization(current, unit.traceAgc); + break; + case ProcStepType::StepHilbertTransform: + current = hilbertTransform(current, unit.hilbert); + break; + case ProcStepType::StepGlobalGain: + current = applyGain(current, unit.gain); + break; + case ProcStepType::StepMigration: + current = kirchhoffMigration(current, unit.migration); + break; + default: + // 未知步骤类型跳过 + break; + } + + qDebug() << "[RadarProcessor] Step" << (i + 1) << "finished:" << stepName(unit.type); + + if (context && context->reportProgress) + context->reportProgress(i + 1, stepCount, stepName(unit.type)); + + if (canceled()) + break; + } + return current; +} \ No newline at end of file diff --git a/external/gpr3dviewer/RadarProcessor.h b/external/gpr3dviewer/RadarProcessor.h new file mode 100644 index 0000000..6983f5a --- /dev/null +++ b/external/gpr3dviewer/RadarProcessor.h @@ -0,0 +1,346 @@ +/* + * RadarProcessor.h + * 探地雷达GPR数据预处理算法核心类头文件 + * 功能：封装全套雷达道数据校正、滤波、增益、均衡、希尔伯特包络算法 + * 架构：模块化单步骤参数结构体 + 流水线顺序调度执行 + * 对接：MainWindow UI界面参数绑定、GPRDataModel三维道集数据模型 + * 适配：RAD/RD3格式多通道、多测线原始雷达数据 + */ +#ifndef RADARPROCESSOR_H +#define RADARPROCESSOR_H + +#include "GPRDataModel.h" +#include +#include +#include + +class RadarProcessor +{ +public: + /** + * @brief 所有可执行预处理步骤枚举标识 + * 流水线严格按照枚举顺序业务逻辑排布 + */ + enum class ProcStepType + { + StepTimeZeroCut, // 0：零时校正（切除直达波零点） + StepZeroDrift, // 1：道基线零漂消除 + StepRemoveDC, // 2：简易整道直流偏移扣除 + StepBackgroundRemove, // 3：剖面背景均值压制 + StepSphericalTVG, // 4：球面扩散+介质吸收TVG深度增益 + StepBandpassFilter, // 5：FIR带通滤波剔除高低频噪声 + StepProfileSmooth, // 6：时空剖面平滑降噪 + StepInnerAGC, // 7：单道内自适应振幅均衡AGC + StepTraceToTraceAGC, // 8：道与道之间整体振幅均衡 + StepHilbertTransform, // 9：希尔伯特变换（包络/正交分量） + StepGlobalGain, // 10：全局统一倍率增益放大 + StepMigration // 11：Kirchhoff偏移 + }; + + //========================================================================= + // 1. 零时切除参数结构体 StepTimeZeroCut + // 作用：定位直达波起跳零点，裁掉零点前无效采样点，统一道起始位置 + //========================================================================= + struct TimeZeroCutParams + { + bool autoDetect = true; // true=自动识别零点；false=手动固定裁掉前N个采样 + int cutSamples = 30; // 手动模式裁切采样点数 + int frontSearchWindow = 180; // 零点搜索窗口：只在道前180个采样内找起跳点 + double noiseSigmaMultiple = 3.0; // 起跳判定阈值：噪声标准差倍数，大于判定为有效直达波 + bool enable = false; // 流水线开关：是否启用本步骤 + }; + + /** + * @brief 零漂消除模式枚举 + */ + enum class ZeroDriftMode + { + DC, // 模式1：整道均值直流偏移一次性扣除（速度快） + Sliding // 模式2：滑动窗口逐段扣除基线（适合长时漂移严重数据） + }; + struct ZeroDriftParams + { + ZeroDriftMode mode = ZeroDriftMode::Sliding; + int slidingWindowSamples = 100; // 滑动窗口采样点数，仅Sliding模式生效 + bool enable = false; + }; + + //========================================================================= + // TVG球面深度增益 StepSphericalTVG + // 物理原理：电磁波随传播距离球面扩散衰减+介质吸收衰减 + //========================================================================= + struct SphericalTvgParams { + bool enableSpherical = true; // 开启球面扩散补偿 + bool enableAbsorption = true; // 开启介质吸收衰减补偿 + double velocityMPerNs = 0.12; // 雷达波速(m/ns)，从GPRDataHeader自动读取 + double referenceDepthM = 0.01; // 参考基准深度（增益归一化基准） + double exponent = 1.5; // 扩散指数：空气1.0、土体/混凝土常用1.0~2.0 + double absorptionBeta = 1.0; // 吸收系数 β (dB/m)，干土小、湿土大 + double maxGain = 30.0; // 最大增益上限，防止深层噪声过度放大 + bool enable = false; + }; + + //========================================================================= + // 带通滤波 StepBandpassFilter + // 滤除天线主频外低频漂移、高频仪器白噪声 + //========================================================================= + struct BandpassParams + { + bool autoFreq = true; // true自动根据天线主频计算通带；false手动填高低频 + double lowFreqHz = 1000; // 通带下限(Hz) + double highFreqHz = 100; // 通带上限(Hz) + double antennaFreqMHz = 200.0; // 天线中心频率(MHz)，自动模式读取头文件 + bool enable = false; + }; + + //========================================================================= + // 剖面二维平滑 StepProfileSmooth + // 垂直=深度方向平滑；水平=测线道方向平滑 + //========================================================================= + struct ProfileSmoothParams + { + int smoothWindow = 2; // 半窗口大小，实际总窗口=2*window+1 + bool verticalSmooth = true; // 深度方向平滑开关 + bool horizontalSmooth = true; // 测线道方向平滑开关 + bool enable = false; + }; + + //========================================================================= + // 背景去除 StepBackgroundRemove + // 压制整剖面静态杂波、地面耦合固定反射、钢筋/管线持续干扰 + //========================================================================= + enum class BackgroundMode + { + MeanAverage, // 均值法：多道平均作为背景相减，通用稳定 + SingularityFilter // 奇异值滤波：保留强异常，压制平缓背景（空洞/管线优选） + }; + struct BackgroundParams + { + BackgroundMode mode = BackgroundMode::MeanAverage; + int averageTraceCount = 301; // 参与平均背景的同通道道数量 + double singularityThreshold = 1.8; // 奇异判定阈值，越大越不容易抹掉小异常 + bool enable = false; + }; + + //========================================================================= + // 道内AGC StepInnerAGC + // 同一道内深浅振幅均衡，消除天然深度衰减（TVG补充均衡） + //========================================================================= + struct TraceInnerAgcParams { + int windowSamples = 120; // 道内滑动RMS统计窗口 + double targetRms = 0.0; // 目标RMS值；0=自适应均衡 + double maxGain = 6.0; // 单窗口最大增益限制 + double epsilon = 1.0; // 防除零极小值 + bool enable = false; + }; + + //========================================================================= + // 道间AGC StepTraceToTraceAGC + // 整条测线所有道之间振幅拉平，消除收发耦合波动、行走速度不均 + //========================================================================= + enum class TraceToTraceMode + { + Global, // 全局：全部道统一均衡基准 + Local // 局部：滑动窗口相邻道均衡（起伏大路面优选） + }; + struct TraceToTraceAgcParams { + TraceToTraceMode mode = TraceToTraceMode::Local; + int horizontalWindowTraces = 31; // 局部模式滑动道窗口 + double targetRms = 0.0; + double maxGain = 4.0; + double epsilon = 1.0; + bool enable = false; + }; + + //========================================================================= + // 希尔伯特变换 StepHilbertTransform + // 输出振幅包络（成像常用）或正交虚部（相位分析） + //========================================================================= + struct HilbertParams + { + bool computeEnvelope = true; // true=振幅包络；false=正交相位分量 + bool enable = false; + }; + + // 全局倍率增益 + struct GlobalGainParams + { + double factor = 1.0; // 放大倍数，0.1~10区间常用 + bool enable = false; + }; + + // 简易去直流（轻量零漂） + struct DcShiftParams + { + bool enable = false; + }; + + //========================================================================= + // Kirchhoff偏移参数 + //========================================================================= + struct MigrationParams + { + int sumWidth = 64; // 两侧各求和的迹线数量 + double velocityMPerNs = 0.12; // 雷达波速(m/ns) + bool enable = false; + }; + + //========================================================================= + // 流水线单步单元：存储步骤类型 + 全套独立参数 + //========================================================================= + struct ProcStepUnit + { + ProcStepType type; + + // 全部算法参数容器，运行时仅type对应的结构体生效 + TimeZeroCutParams zeroCut; + SphericalTvgParams tvg; + BandpassParams band; + ProfileSmoothParams smooth; + BackgroundParams bg; + ZeroDriftParams zeroDrift; + TraceInnerAgcParams innerAgc; + TraceToTraceAgcParams traceAgc; + HilbertParams hilbert; + GlobalGainParams gain; + DcShiftParams dc; + MigrationParams migration; + }; + + //========================================================================= + // 完整预处理流水线：有序步骤队列 + 默认标准流程 + //========================================================================= + struct ProcPipeline + { + QVector steps; + + /** + * @brief 加载工业通用标准处理流水线（道路检测默认） + * 顺序：零时→零漂→背景→TVG增益→带通→平滑→内AGC→道间AGC→包络 + */ + void setDefaultFlow() + { + steps.clear(); + + // 1. 零时校正 + ProcStepUnit s0; + s0.type = ProcStepType::StepTimeZeroCut; + s0.zeroCut.autoDetect = true; + s0.zeroCut.cutSamples = 30; + s0.zeroCut.frontSearchWindow = 180; + s0.zeroCut.noiseSigmaMultiple = 3.0; + s0.zeroCut.enable = true; + steps.append(s0); + + // 2. 零漂去除 + ProcStepUnit s1; + s1.type = ProcStepType::StepZeroDrift; + s1.zeroDrift.mode = ZeroDriftMode::Sliding; + s1.zeroDrift.slidingWindowSamples = 100; + s1.zeroDrift.enable = true; + steps.append(s1); + + // 3. 背景压制 + ProcStepUnit s2; + s2.type = ProcStepType::StepBackgroundRemove; + s2.bg.mode = BackgroundMode::MeanAverage; + s2.bg.averageTraceCount = 301; + s2.bg.enable = true; + steps.append(s2); + + // 4. TVG深度增益 + ProcStepUnit s3; + s3.type = ProcStepType::StepSphericalTVG; + s3.tvg.enableSpherical = true; + s3.tvg.enableAbsorption = true; + s3.tvg.exponent = 1.0; + s3.tvg.maxGain = 20.0; + s3.tvg.absorptionBeta = 0.002; + s3.tvg.enable = true; + steps.append(s3); + + // 5. 带通滤波 + ProcStepUnit s4; + s4.type = ProcStepType::StepBandpassFilter; + s4.band.autoFreq = true; + s4.band.lowFreqHz = 100; + s4.band.highFreqHz = 1500; + s4.band.enable = true; + steps.append(s4); + + // 6. 剖面平滑 + ProcStepUnit s5; + s5.type = ProcStepType::StepProfileSmooth; + s5.smooth.smoothWindow = 2; + s5.smooth.verticalSmooth = true; + s5.smooth.horizontalSmooth = true; + s5.smooth.enable = true; + steps.append(s5); + + // 7. 道内AGC + ProcStepUnit s6; + s6.type = ProcStepType::StepInnerAGC; + s6.innerAgc.windowSamples = 120; + s6.innerAgc.maxGain = 6.0; + s6.innerAgc.enable = true; + steps.append(s6); + + // 8. 道间均衡AGC + ProcStepUnit s7; + s7.type = ProcStepType::StepTraceToTraceAGC; + s7.traceAgc.mode = TraceToTraceMode::Local; + s7.traceAgc.horizontalWindowTraces = 31; + s7.traceAgc.maxGain = 4.0; + s7.traceAgc.enable = true; + steps.append(s7); + } + }; + +public: + struct PipelineRuntimeContext { + std::function isCanceled; + std::function reportProgress; + }; + + RadarProcessor(); + + //==================== 单算法静态处理函数 ==================== + /// 简易去除整道直流偏移 + static GPRDataModel removeDCShift(const GPRDataModel &input, const DcShiftParams& param); + /// 零漂基线校正 + static GPRDataModel removeZeroDrift(const GPRDataModel &input, const ZeroDriftParams& param); + /// 全局统一倍率增益 + static GPRDataModel applyGain(const GPRDataModel &input, const GlobalGainParams& param); + /// TVG球面+吸收深度增益补偿 + static GPRDataModel sphericalTvg(const GPRDataModel &input, const SphericalTvgParams ¶ms); + /// 二维剖面均值平滑 + static GPRDataModel profileSmooth(const GPRDataModel &input, const ProfileSmoothParams ¶ms); + /// 道内滑动AGC均衡 + static GPRDataModel traceInnerAgc(const GPRDataModel &input, const TraceInnerAgcParams ¶ms); + /// 道与道之间振幅均衡 + static GPRDataModel traceToTraceEqualization(const GPRDataModel &input, const TraceToTraceAgcParams ¶ms); + /// 希尔伯特变换包络/相位 + static GPRDataModel hilbertTransform(const GPRDataModel &input, const HilbertParams ¶ms); + + /// 零点自动/手动裁切算法实现 + static GPRDataModel cutTimeZero(const GPRDataModel& input, const TimeZeroCutParams& param); + /// FIR带通滤波实现 + static GPRDataModel bandpassFilter(const GPRDataModel &input, const BandpassParams ¶ms); + /// 背景杂波去除实现（均值/奇异滤波双模式） + static GPRDataModel backgroundRemoval(const GPRDataModel &input, const BackgroundParams ¶ms); + /// Kirchhoff偏移（x-t域双曲线求和） + static GPRDataModel kirchhoffMigration(const GPRDataModel &input, const MigrationParams ¶ms); + + /** + * @brief 流水线总执行入口 + * @param rawData 原始未处理GPR三维道集数据 + * @param pipeline 有序步骤流水线配置 + * @return 逐步骤处理完成后的成品数据 + */ + static GPRDataModel runPipeline(const GPRDataModel& rawData, const ProcPipeline& pipeline); + static GPRDataModel runPipeline(const GPRDataModel& rawData, + const ProcPipeline& pipeline, + const PipelineRuntimeContext *context); +}; + +#endif // RADARPROCESSOR_H \ No newline at end of file diff --git a/external/gpr3dviewer/RadarTypes.h b/external/gpr3dviewer/RadarTypes.h new file mode 100644 index 0000000..a9e069d --- /dev/null +++ b/external/gpr3dviewer/RadarTypes.h @@ -0,0 +1,9 @@ +#ifndef RADARTYPES_H +#define RADARTYPES_H + +enum class RadarType { + Mala_Mira, // .rad + .rd3 + Impulse // .iprh + .iprb +}; + +#endif // RADARTYPES_H diff --git a/external/gpr3dviewer/Rd3Parser.cpp b/external/gpr3dviewer/Rd3Parser.cpp new file mode 100644 index 0000000..3990795 --- /dev/null +++ b/external/gpr3dviewer/Rd3Parser.cpp @@ -0,0 +1,300 @@ +#include "Rd3Parser.h" +#include "PerformanceLogger.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/** + * @brief 加载整套Mala Mira系列雷达RD3数据（rad头文件 + rd3纯二进制波形） + * @param radFilePath .rad文本配置头文件路径 + * @param model 输出GPR全局数据模型 + * @return true 加载解析成功；false 失败 + * @note 存储结构：.rad存储全部仪器/测线参数；.rd3无文件头，从头到尾全是short16振幅采样 + */ +bool Rd3Parser::loadFromRad(const QString &radFilePath, GPRDataModel &model) { + SCOPED_PERF_TIMER("Parser", "Rd3Parser::loadFromRad"); + + // 清空旧数据与头部信息 + model.clear(); + model.header = GPRDataModel::Header{}; + + // 第一步：解析rad文本头，填充全部采集参数 + if (!parseRadHeader(radFilePath, model.header)) { + qDebug() << "Error: Failed to parse .rad header file:" << radFilePath; + return false; + } + + // 自动匹配同目录同名二进制数据文件 .rd3 + QFileInfo radInfo(radFilePath); + QString binaryPath = radInfo.absolutePath() + "/" + radInfo.completeBaseName() + ".rd3"; + + if (!QFile::exists(binaryPath)) { + qDebug() << "Error: Matching binary file not found at:" << binaryPath; + return false; + } + + // 推算后仍无有效道数则失败 + if (model.header.numTraces <= 0) { + qDebug() << "Error: Unable to determine numTraces from header"; + return false; + } + + // 第二步：读取无头部的纯二进制波形 + return loadRd3Binary(binaryPath, model); +} + +/** + * @brief 解析RAD文本头文件，提取雷达采集关键参数存入Header结构体 + * @param radFilePath rad文本文件路径 + * @param header 待填充头部参数结构体 + * @return 解析成功返回true + */ +bool Rd3Parser::parseRadHeader(const QString &radFilePath, GPRDataModel::Header &header) { + SCOPED_PERF_TIMER("Parser", "Rd3Parser::parseRadHeader"); + + QFile file(radFilePath); + // 只读文本模式打开 + if (!file.open(QIODevice::ReadOnly | QIODevice::Text)) { + qDebug() << "Error: Cannot open rad file for read"; + return false; + } + + QTextStream in(&file); + + // 逐行读取键值对 KEY: VALUE 或 KEY=VALUE + while (!in.atEnd()) { + QString line = in.readLine().trimmed(); + if (line.isEmpty()) continue; + + int sepPos = line.indexOf(':'); + if (sepPos == -1) sepPos = line.indexOf('='); + if (sepPos == -1) continue; + + QString key = line.left(sepPos).trimmed(); + QString value = line.mid(sepPos + 1).trimmed(); + + // 原始字符串参数全存入map备用 + header.rawParams[key] = value; + + // 识别核心业务参数并类型转换赋值 + if (key == "SAMPLES") { + header.samplesPerTrace = extractInt(value); + } else if (key == "LAST TRACE") { + header.numTraces = extractInt(value); + } else if (key == "TIMEWINDOW") { + header.timeWindowNs = extractDouble(value); + } else if (key == "DISTANCE INTERVAL") { + header.distanceInc = extractDouble(value); + } else if (key == "ANTENNAS") { + header.antennaFreq = extractDouble(value); + } else if (key == "ANTENNAS") { + header.antennaType = value; + } else if (key == "DATE") { + header.date = value; + } else if (key == "TIME") { + header.timeStr = value; + } else if (key == "NUMBER_OF_CH") { + header.numberOfChannels = extractInt(value); + } else if (key == "CH_X_OFFSETS") { + // 多通道天线水平X偏移量 + QStringList offsets = value.split(' ', Qt::SkipEmptyParts); + for (const QString &offset : offsets) { + header.chXOffsets.append(offset.toFloat()); + } + } else if (key == "CH_Y_OFFSETS") { + // 多通道天线行进Y偏移量 + QStringList offsets = value.split(' ', Qt::SkipEmptyParts); + for (const QString &offset : offsets) { + header.chYOffsets.append(offset.toFloat()); + } + } else if (key == "UNITS") { + header.units = value; + } else if (key == "START POSITION") { + header.startPosition = extractDouble(value); + } else if (key == "STOP POSITION") { + header.stopPosition = extractDouble(value); + } else if (key == "TIME INTERVAL") { + header.timeIntervalNs = extractDouble(value); + } + } + + file.close(); + + // 基础合法性校验：单道采样数必须大于0；总道数允许从二进制文件推算 + if (header.samplesPerTrace <= 0) { + qDebug() << "Error: Invalid SAMPLES value in rad file"; + return false; + } + + // 设置地层电磁波默认速度 0.1 m/ns（空气近似值，后续界面可手动修改） + header.waveVelocity = 0.1; + + qDebug() << "==== RAD Header Parse Complete ====" + << "\nSamples per trace:" << header.samplesPerTrace + << "\nTotal traces:" << header.numTraces + << "\nTime window(ns):" << header.timeWindowNs + << "\nChannel count:" << header.numberOfChannels + << "\nX offset array size:" << header.chXOffsets.size() + << "\nY offset array size:" << header.chYOffsets.size(); + + return true; +} + +/** + * @brief 读取纯RD3二进制数据（无任何文件头，全连续short16振幅） + * @param rd3FilePath rd3波形文件路径 + * @param model 绑定头部参数，填充完整traces波形数组 + * @return 读取成功true + * @detail 存储格式： + * 道0采样0,道0采样1...道0采样N | 道1采样0...道1采样N | ...循环所有通道所有道 + * 数据类型：小端序 short 16位有符号整型 + */ +bool Rd3Parser::loadRd3Binary(const QString &rd3FilePath, GPRDataModel &model) { + SCOPED_PERF_TIMER("Parser", "Rd3Parser::loadRd3Binary"); + + QFile file(rd3FilePath); + if (!file.open(QIODevice::ReadOnly)) { + qDebug() << "Error: Open rd3 binary failed:" << rd3FilePath; + return false; + } + + const int samplesPerTrace = model.header.samplesPerTrace; + const int totalTraceCount = model.header.numTraces; + + // ========== 关键修正：rd3无头部偏移，直接从文件0位置开始读波形 ========== + const qint64 dataStartOffset = 0; + file.seek(dataStartOffset); + + // 理论完整字节大小校验 + const qint64 singleTraceByteSize = samplesPerTrace * sizeof(short); + const qint64 fullExpectedBytes = totalTraceCount * singleTraceByteSize; + const qint64 realFileBytes = file.size(); + + if (realFileBytes < fullExpectedBytes) { + qDebug() << "Warning: RD3 file size smaller than theoretical data size!" + << "Expected:" << fullExpectedBytes << "Actual:" << realFileBytes; + } + + // 预分配内存容器，减少动态扩容开销 + model.traces.reserve(totalTraceCount); + + // 单道读取缓冲区 + QByteArray traceBuffer; + traceBuffer.resize(singleTraceByteSize); + + // 通道数量兜底 + const int channelCnt = model.header.numberOfChannels > 0 ? model.header.numberOfChannels : 1; + model.channels = channelCnt; + // 每通道独立道数量 + model.tracesPerChannel = totalTraceCount / channelCnt; + + // 计算整条测线总行进距离 + if (model.header.distanceInc > 1e-6) { + model.totalDistance = static_cast(model.tracesPerChannel * model.header.distanceInc); + } else { + model.totalDistance = static_cast(model.header.stopPosition - model.header.startPosition); + } + + // 循环逐道读取二进制波形 + for (int traceGlobalIdx = 0; traceGlobalIdx < totalTraceCount; ++traceGlobalIdx) { + + if (file.atEnd()) { + qDebug() << "Warning: File ended early at global trace index" << traceGlobalIdx; + break; + } + + // 读取一整道所有采样字节 + qint64 readBytes = file.read(traceBuffer.data(), traceBuffer.size()); + if (readBytes != traceBuffer.size()) { + qDebug() << "Warning: Trace" << traceGlobalIdx << "incomplete byte read"; + break; + } + + RadarTrace oneTrace; + oneTrace.amplitudes.resize(samplesPerTrace); + // 二进制指针强转short（Windows平台原生小端序，雷达标准字节序） + const short* rawShortBuf = reinterpret_cast(traceBuffer.constData()); + + // 填充单道振幅值 + for (int s = 0; s < samplesPerTrace; s++) { + oneTrace.amplitudes[s] = rawShortBuf[s]; + } + + // 分配所属通道、通道内道序号 + int chNo = traceGlobalIdx % channelCnt; + int traceInChIdx = traceGlobalIdx / channelCnt; + oneTrace.channelNumber = chNo; + + // 计算空间三维坐标 X(天线横向偏移) Y(行进里程) Z=0地面平面 + float xOff = 0.0f; + float yOff = 0.0f; + if (chNo < model.header.chXOffsets.size()) xOff = model.header.chXOffsets[chNo]; + if (chNo < model.header.chYOffsets.size()) yOff = model.header.chYOffsets[chNo]; + + float lineDist = static_cast(model.header.startPosition + traceInChIdx * model.header.distanceInc); + oneTrace.position = QVector3D(xOff, lineDist - yOff, 0.0f); + + model.traces.append(std::move(oneTrace)); + } + + file.close(); + + if (model.traces.isEmpty()) { + qDebug() << "Error: No valid traces loaded from rd3"; + return false; + } + + qDebug() << "RD3 Binary Load OK, total valid traces:" << model.traces.size(); + return true; +} + +/** + * @brief 安全字符串转double数值，优先提取字符串开头的数字，转换失败返回0 + * @param value 原始字符串（如：200 MHz shielded） + * @return 浮点结果 + */ +double Rd3Parser::extractDouble(const QString& value) { + bool ok = false; + double res = 0.0; + + // 1. 遍历字符串，截取开头连续的数字/小数点部分 + int numLength = 0; + while (numLength < value.length()) { + QChar ch = value.at(numLength); + // 只保留 数字 和 小数点（支持小数，如 200.5 MHz） + if (ch.isDigit() || ch == '.') { + numLength++; + } + else { + // 遇到非数字/小数点，停止截取 + break; + } + } + + // 2. 截取有效数字字符串并转换 + if (numLength > 0) { + QString numStr = value.left(numLength); + res = numStr.toDouble(&ok); + } + + // 3. 转换成功返回数值，失败返回0.0 + return ok ? res : 0.0; +} + +/** + * @brief 安全字符串转int整型，转换失败返回0 + * @param value 原始rad文件字符串数值 + * @return 整型结果 + */ +int Rd3Parser::extractInt(const QString &value) { + bool ok = false; + int res = value.toInt(&ok); + return ok ? res : 0; +} + diff --git a/external/gpr3dviewer/Rd3Parser.h b/external/gpr3dviewer/Rd3Parser.h new file mode 100644 index 0000000..a21c990 --- /dev/null +++ b/external/gpr3dviewer/Rd3Parser.h @@ -0,0 +1,21 @@ +#ifndef RD3PARSER_H +#define RD3PARSER_H + +#include "GPRDataModel.h" +#include + +class Rd3Parser { +public: + // 主入口：传入 .rad 文件路径，它会自动寻找同名的 .rd3 文件 + static bool loadFromRad(const QString &radFilePath, GPRDataModel &model); + +private: + static bool parseRadHeader(const QString &radFilePath, GPRDataModel::Header &header); + static bool loadRd3Binary(const QString &rd3FilePath, GPRDataModel &model); + + // 辅助：从 .rad 的字符串中提取数值 + static double extractDouble(const QString &value); + static int extractInt(const QString &value); +}; + +#endif // RD3PARSER_H \ No newline at end of file diff --git a/external/gpr3dviewer/SurveyGeometry.h b/external/gpr3dviewer/SurveyGeometry.h new file mode 100644 index 0000000..deeb6f1 --- /dev/null +++ b/external/gpr3dviewer/SurveyGeometry.h @@ -0,0 +1,127 @@ +#ifndef SURVEYGEOMETRY_H +#define SURVEYGEOMETRY_H + +#include +#include +#include + +struct SurveyGeometry { + double rtkOffsetX = 0.0; // RTK相对天线中心的设备x轴偏移（m），头文件无对应字段，保留给用户输入 + double rtkOffsetY = 0.0; // RTK相对天线中心的设备y轴偏移（m），来自CH_Y_OFFSETS + double ch1XRel = 0.0; // 第1通道相对天线中心的x坐标（m），由CH_X_OFFSETS首尾值推导 + QVector channelXRel; // 各通道相对天线中心的x坐标（m），来自CH_X_OFFSETS归一到天线中心 + int channelCount = 1; // 通道总数，来自NUMBER_OF_CH + int cgcsZone = 0; // CGCS2000 3度带带号，0表示自动检测 + double centralMeridianDeg = 0.0; // 中央经线（度），自动推导 + + bool operator==(const SurveyGeometry &other) const { + return rtkOffsetX == other.rtkOffsetX && + rtkOffsetY == other.rtkOffsetY && + ch1XRel == other.ch1XRel && + channelXRel == other.channelXRel && + channelCount == other.channelCount && + cgcsZone == other.cgcsZone && + centralMeridianDeg == other.centralMeridianDeg; + } + + void applyHeaderOffsets(int numberOfChannels, + const QVector &chXOffsets, + const QVector &chYOffsets) + { + const int offsetChannelCount = qMax(chXOffsets.size(), chYOffsets.size()); + channelCount = qMax(1, numberOfChannels > 0 ? numberOfChannels : offsetChannelCount); + rtkOffsetY = chYOffsets.isEmpty() ? 0.0 : static_cast(chYOffsets.first()); + + channelXRel.clear(); + channelXRel.reserve(channelCount); + + if (chXOffsets.size() >= 2) { + const double centerX = (static_cast(chXOffsets.first()) + chXOffsets.last()) / 2.0; + for (int i = 0; i < chXOffsets.size() && i < channelCount; ++i) { + channelXRel.append(static_cast(chXOffsets[i]) - centerX); + } + ch1XRel = channelXRel.isEmpty() ? 0.0 : channelXRel.first(); + } else { + ch1XRel = 0.0; + } + + if (channelXRel.size() < channelCount) { + const int existing = channelXRel.size(); + channelXRel.resize(channelCount); + if (existing == 0) { + for (int c = 0; c < channelCount; ++c) { + channelXRel[c] = 0.0; + } + } else if (channelCount > 1) { + const double lastXRel = -ch1XRel; + for (int c = existing; c < channelCount; ++c) { + channelXRel[c] = ch1XRel + (lastXRel - ch1XRel) * c / (channelCount - 1.0); + } + } + } + } + + static SurveyGeometry fromHeaderOffsets(int numberOfChannels, + const QVector &chXOffsets, + const QVector &chYOffsets) + { + SurveyGeometry g; + g.applyHeaderOffsets(numberOfChannels, chXOffsets, chYOffsets); + return g; + } + + QJsonObject toJson() const { + QJsonObject obj; + obj["rtkOffsetX"] = rtkOffsetX; + obj["rtkOffsetY"] = rtkOffsetY; + obj["ch1XRel"] = ch1XRel; + QJsonArray channelXRelArray; + for (double x : channelXRel) { + channelXRelArray.append(x); + } + obj["channelXRel"] = channelXRelArray; + obj["channelCount"] = channelCount; + obj["cgcsZone"] = cgcsZone; + obj["centralMeridianDeg"] = centralMeridianDeg; + return obj; + } + + static SurveyGeometry fromJson(const QJsonObject &obj) { + SurveyGeometry g; + g.rtkOffsetX = obj.value("rtkOffsetX").toDouble(0.0); + g.rtkOffsetY = obj.value("rtkOffsetY").toDouble(0.0); + g.ch1XRel = obj.value("ch1XRel").toDouble(0.0); + g.channelCount = qMax(1, obj.value("channelCount").toInt(1)); + g.cgcsZone = obj.value("cgcsZone").toInt(0); + g.centralMeridianDeg = obj.value("centralMeridianDeg").toDouble(0.0); + + const QJsonArray channelXRelArray = obj.value("channelXRel").toArray(); + for (const QJsonValue &value : channelXRelArray) { + g.channelXRel.append(value.toDouble(0.0)); + } + + if (g.channelXRel.isEmpty() && obj.contains("ch16XRel")) { + const double legacyLastXRel = obj.value("ch16XRel").toDouble(-g.ch1XRel); + g.channelXRel.resize(g.channelCount); + for (int c = 0; c < g.channelCount; ++c) { + g.channelXRel[c] = (g.channelCount > 1) + ? g.ch1XRel + (legacyLastXRel - g.ch1XRel) * c / (g.channelCount - 1.0) + : g.ch1XRel; + } + } else if (g.channelXRel.isEmpty()) { + g.channelXRel.resize(g.channelCount); + const double lastXRel = -g.ch1XRel; + for (int c = 0; c < g.channelCount; ++c) { + g.channelXRel[c] = (g.channelCount > 1) + ? g.ch1XRel + (lastXRel - g.ch1XRel) * c / (g.channelCount - 1.0) + : g.ch1XRel; + } + } else if (g.channelXRel.size() != g.channelCount) { + g.channelCount = qMax(1, g.channelXRel.size()); + } + + return g; + } +}; + +#endif // SURVEYGEOMETRY_H diff --git a/external/gpr3dviewer/third_party/kissfft/_kiss_fft_guts.h b/external/gpr3dviewer/third_party/kissfft/_kiss_fft_guts.h new file mode 100644 index 0000000..f4b6832 --- /dev/null +++ b/external/gpr3dviewer/third_party/kissfft/_kiss_fft_guts.h @@ -0,0 +1,189 @@ +/* + * Copyright (c) 2003-2010, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +/* kiss_fft.h + defines kiss_fft_scalar as either short or a float type + and defines + typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */ + +#ifndef _kiss_fft_guts_h +#define _kiss_fft_guts_h + +#include "kiss_fft.h" +#include "kiss_fft_log.h" +#include + +#define MAXFACTORS 32 +/* e.g. an fft of length 128 has 4 factors + as far as kissfft is concerned + 4*4*4*2 + */ + +struct kiss_fft_state{ + int nfft; + int inverse; + int factors[2*MAXFACTORS]; + kiss_fft_cpx twiddles[1]; +}; + +/* + Explanation of macros dealing with complex math: + + C_MUL(m,a,b) : m = a*b + C_FIXDIV( c , div ) : if a fixed point impl., c /= div. noop otherwise + C_SUB( res, a,b) : res = a - b + C_SUBFROM( res , a) : res -= a + C_ADDTO( res , a) : res += a + * */ +#ifdef FIXED_POINT +#include +#if (FIXED_POINT==32) +# define FRACBITS 31 +# define SAMPPROD int64_t +#define SAMP_MAX INT32_MAX +#define SAMP_MIN INT32_MIN +#else +# define FRACBITS 15 +# define SAMPPROD int32_t +#define SAMP_MAX INT16_MAX +#define SAMP_MIN INT16_MIN +#endif + +#if defined(CHECK_OVERFLOW) +# define CHECK_OVERFLOW_OP(a,op,b) \ + if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \ + KISS_FFT_WARNING("overflow (%d " #op" %d) = %ld", (a),(b),(SAMPPROD)(a) op (SAMPPROD)(b)); } +#endif + + +# define smul(a,b) ( (SAMPPROD)(a)*(b) ) +# define sround( x ) (kiss_fft_scalar)( ( (x) + (1<<(FRACBITS-1)) ) >> FRACBITS ) + +# define S_MUL(a,b) sround( smul(a,b) ) + +# define C_MUL(m,a,b) \ + do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \ + (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0) + +# define DIVSCALAR(x,k) \ + (x) = sround( smul( x, SAMP_MAX/k ) ) + +# define C_FIXDIV(c,div) \ + do { DIVSCALAR( (c).r , div); \ + DIVSCALAR( (c).i , div); }while (0) + +# define C_MULBYSCALAR( c, s ) \ + do{ (c).r = sround( smul( (c).r , s ) ) ;\ + (c).i = sround( smul( (c).i , s ) ) ; }while(0) + +#else /* not FIXED_POINT*/ + +# define S_MUL(a,b) ( (a)*(b) ) +#define C_MUL(m,a,b) \ + do{ (m).r = (a).r*(b).r - (a).i*(b).i;\ + (m).i = (a).r*(b).i + (a).i*(b).r; }while(0) +# define C_FIXDIV(c,div) /* NOOP */ +# define C_MULBYSCALAR( c, s ) \ + do{ (c).r *= (s);\ + (c).i *= (s); }while(0) +#endif + +#ifndef CHECK_OVERFLOW_OP +# define CHECK_OVERFLOW_OP(a,op,b) /* noop */ +#endif + +#define C_ADD( res, a,b)\ + do { \ + CHECK_OVERFLOW_OP((a).r,+,(b).r)\ + CHECK_OVERFLOW_OP((a).i,+,(b).i)\ + (res).r=(a).r+(b).r; (res).i=(a).i+(b).i; \ + }while(0) +#define C_SUB( res, a,b)\ + do { \ + CHECK_OVERFLOW_OP((a).r,-,(b).r)\ + CHECK_OVERFLOW_OP((a).i,-,(b).i)\ + (res).r=(a).r-(b).r; (res).i=(a).i-(b).i; \ + }while(0) +#define C_ADDTO( res , a)\ + do { \ + CHECK_OVERFLOW_OP((res).r,+,(a).r)\ + CHECK_OVERFLOW_OP((res).i,+,(a).i)\ + (res).r += (a).r; (res).i += (a).i;\ + }while(0) + +#define C_SUBFROM( res , a)\ + do {\ + CHECK_OVERFLOW_OP((res).r,-,(a).r)\ + CHECK_OVERFLOW_OP((res).i,-,(a).i)\ + (res).r -= (a).r; (res).i -= (a).i; \ + }while(0) + + +#ifdef FIXED_POINT +# define KISS_FFT_COS(phase) floor(.5+SAMP_MAX * cos (phase)) +# define KISS_FFT_SIN(phase) floor(.5+SAMP_MAX * sin (phase)) +# define HALF_OF(x) ((x)>>1) +#elif defined(USE_SIMD) +#if defined(HAVE_LASX) +#define KISS_FFT_COS(phase) ({ \ + float __cos_val = cosf(phase); \ + (__m256)(__lasx_xvldrepl_w(&__cos_val, 0)); \ +}) +#define KISS_FFT_SIN(phase) ({ \ + float __sin_val = sinf(phase); \ + (__m256)(__lasx_xvldrepl_w(&__sin_val, 0)); \ +}) +#define HALF_OF(x) ((x) * (__m256)(__lasx_xvreplgr2vr_w(0x3F000000))) // 0.5f +#elif defined(HAVE_LSX) +#define KISS_FFT_COS(phase) ({ \ + float __cos_val = cosf(phase); \ + (__m128)(__lsx_vldrepl_w(&__cos_val, 0)); \ +}) +#define KISS_FFT_SIN(phase) ({ \ + float __sin_val = sinf(phase); \ + (__m128)(__lsx_vldrepl_w(&__sin_val, 0)); \ +}) +#define HALF_OF(x) ((x) * (__m128)(__lsx_vreplgr2vr_w(0x3F000000))) // 0.5f +#else +# define KISS_FFT_COS(phase) _mm_set1_ps( cos(phase) ) +# define KISS_FFT_SIN(phase) _mm_set1_ps( sin(phase) ) +# define HALF_OF(x) ((x)*_mm_set1_ps(.5)) +#endif +#else +# define KISS_FFT_COS(phase) (kiss_fft_scalar) cos(phase) +# define KISS_FFT_SIN(phase) (kiss_fft_scalar) sin(phase) +# define HALF_OF(x) ((x)*((kiss_fft_scalar).5)) +#endif + +#define kf_cexp(x,phase) \ + do{ \ + (x)->r = KISS_FFT_COS(phase);\ + (x)->i = KISS_FFT_SIN(phase);\ + }while(0) + + +/* a debugging function */ +#define pcpx(c)\ + KISS_FFT_DEBUG("%g + %gi\n",(double)((c)->r),(double)((c)->i)) + + +#ifdef KISS_FFT_USE_ALLOCA +// define this to allow use of alloca instead of malloc for temporary buffers +// Temporary buffers are used in two case: +// 1. FFT sizes that have "bad" factors. i.e. not 2,3 and 5 +// 2. "in-place" FFTs. Notice the quotes, since kissfft does not really do an in-place transform. +#include +#define KISS_FFT_TMP_ALLOC(nbytes) alloca(nbytes) +#define KISS_FFT_TMP_FREE(ptr) +#else +#define KISS_FFT_TMP_ALLOC(nbytes) KISS_FFT_MALLOC(nbytes) +#define KISS_FFT_TMP_FREE(ptr) KISS_FFT_FREE(ptr) +#endif + +#endif /* _kiss_fft_guts_h */ + diff --git a/external/gpr3dviewer/third_party/kissfft/kiss_fft.c b/external/gpr3dviewer/third_party/kissfft/kiss_fft.c new file mode 100644 index 0000000..aba63e0 --- /dev/null +++ b/external/gpr3dviewer/third_party/kissfft/kiss_fft.c @@ -0,0 +1,424 @@ +/* + * Copyright (c) 2003-2010, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#include +#include "_kiss_fft_guts.h" +/* The guts header contains all the multiplication and addition macros that are defined for + fixed or floating point complex numbers. It also delares the kf_ internal functions. + */ + +static void kf_bfly2( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + int m + ) +{ + kiss_fft_cpx * Fout2; + kiss_fft_cpx * tw1 = st->twiddles; + kiss_fft_cpx t; + Fout2 = Fout + m; + do{ + C_FIXDIV(*Fout,2); C_FIXDIV(*Fout2,2); + + C_MUL (t, *Fout2 , *tw1); + tw1 += fstride; + C_SUB( *Fout2 , *Fout , t ); + C_ADDTO( *Fout , t ); + ++Fout2; + ++Fout; + }while (--m); +} + +static void kf_bfly4( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + const size_t m + ) +{ + kiss_fft_cpx *tw1,*tw2,*tw3; + kiss_fft_cpx scratch[6]; + size_t k=m; + const size_t m2=2*m; + const size_t m3=3*m; + + + tw3 = tw2 = tw1 = st->twiddles; + + do { + C_FIXDIV(*Fout,4); C_FIXDIV(Fout[m],4); C_FIXDIV(Fout[m2],4); C_FIXDIV(Fout[m3],4); + + C_MUL(scratch[0],Fout[m] , *tw1 ); + C_MUL(scratch[1],Fout[m2] , *tw2 ); + C_MUL(scratch[2],Fout[m3] , *tw3 ); + + C_SUB( scratch[5] , *Fout, scratch[1] ); + C_ADDTO(*Fout, scratch[1]); + C_ADD( scratch[3] , scratch[0] , scratch[2] ); + C_SUB( scratch[4] , scratch[0] , scratch[2] ); + C_SUB( Fout[m2], *Fout, scratch[3] ); + tw1 += fstride; + tw2 += fstride*2; + tw3 += fstride*3; + C_ADDTO( *Fout , scratch[3] ); + + if(st->inverse) { + Fout[m].r = scratch[5].r - scratch[4].i; + Fout[m].i = scratch[5].i + scratch[4].r; + Fout[m3].r = scratch[5].r + scratch[4].i; + Fout[m3].i = scratch[5].i - scratch[4].r; + }else{ + Fout[m].r = scratch[5].r + scratch[4].i; + Fout[m].i = scratch[5].i - scratch[4].r; + Fout[m3].r = scratch[5].r - scratch[4].i; + Fout[m3].i = scratch[5].i + scratch[4].r; + } + ++Fout; + }while(--k); +} + +static void kf_bfly3( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + size_t m + ) +{ + size_t k=m; + const size_t m2 = 2*m; + kiss_fft_cpx *tw1,*tw2; + kiss_fft_cpx scratch[5]; + kiss_fft_cpx epi3; + epi3 = st->twiddles[fstride*m]; + + tw1=tw2=st->twiddles; + + do{ + C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3); + + C_MUL(scratch[1],Fout[m] , *tw1); + C_MUL(scratch[2],Fout[m2] , *tw2); + + C_ADD(scratch[3],scratch[1],scratch[2]); + C_SUB(scratch[0],scratch[1],scratch[2]); + tw1 += fstride; + tw2 += fstride*2; + + Fout[m].r = Fout->r - HALF_OF(scratch[3].r); + Fout[m].i = Fout->i - HALF_OF(scratch[3].i); + + C_MULBYSCALAR( scratch[0] , epi3.i ); + + C_ADDTO(*Fout,scratch[3]); + + Fout[m2].r = Fout[m].r + scratch[0].i; + Fout[m2].i = Fout[m].i - scratch[0].r; + + Fout[m].r -= scratch[0].i; + Fout[m].i += scratch[0].r; + + ++Fout; + }while(--k); +} + +static void kf_bfly5( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + int m + ) +{ + kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; + int u; + kiss_fft_cpx scratch[13]; + kiss_fft_cpx * twiddles = st->twiddles; + kiss_fft_cpx *tw; + kiss_fft_cpx ya,yb; + ya = twiddles[fstride*m]; + yb = twiddles[fstride*2*m]; + + Fout0=Fout; + Fout1=Fout0+m; + Fout2=Fout0+2*m; + Fout3=Fout0+3*m; + Fout4=Fout0+4*m; + + tw=st->twiddles; + for ( u=0; ur += scratch[7].r + scratch[8].r; + Fout0->i += scratch[7].i + scratch[8].i; + + scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r); + scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r); + + scratch[6].r = S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i); + scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i); + + C_SUB(*Fout1,scratch[5],scratch[6]); + C_ADD(*Fout4,scratch[5],scratch[6]); + + scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r); + scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r); + scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i); + scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i); + + C_ADD(*Fout2,scratch[11],scratch[12]); + C_SUB(*Fout3,scratch[11],scratch[12]); + + ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; + } +} + +/* perform the butterfly for one stage of a mixed radix FFT */ +static void kf_bfly_generic( + kiss_fft_cpx * Fout, + const size_t fstride, + const kiss_fft_cfg st, + int m, + int p + ) +{ + int u,k,q1,q; + kiss_fft_cpx * twiddles = st->twiddles; + kiss_fft_cpx t; + int Norig = st->nfft; + + kiss_fft_cpx * scratch = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC(sizeof(kiss_fft_cpx)*p); + if (scratch == NULL){ + KISS_FFT_ERROR("Memory allocation failed."); + return; + } + + for ( u=0; u=Norig) twidx-=Norig; + C_MUL(t,scratch[q] , twiddles[twidx] ); + C_ADDTO( Fout[ k ] ,t); + } + k += m; + } + } + KISS_FFT_TMP_FREE(scratch); +} + +static +void kf_work( + kiss_fft_cpx * Fout, + const kiss_fft_cpx * f, + const size_t fstride, + int in_stride, + int * factors, + const kiss_fft_cfg st + ) +{ + kiss_fft_cpx * Fout_beg=Fout; + const int p=*factors++; /* the radix */ + const int m=*factors++; /* stage's fft length/p */ + const kiss_fft_cpx * Fout_end = Fout + p*m; + +#ifdef _OPENMP + // use openmp extensions at the + // top-level (not recursive) + if (fstride==1 && p<=5 && m!=1) + { + int k; + + // execute the p different work units in different threads +# pragma omp parallel for + for (k=0;k floor_sqrt) + p = n; /* no more factors, skip to end */ + } + n /= p; + *facbuf++ = p; + *facbuf++ = n; + } while (n > 1); +} + +/* + * + * User-callable function to allocate all necessary storage space for the fft. + * + * The return value is a contiguous block of memory, allocated with malloc. As such, + * It can be freed with free(), rather than a kiss_fft-specific function. + * */ +kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem ) +{ + KISS_FFT_ALIGN_CHECK(mem) + + kiss_fft_cfg st=NULL; + // check for overflow condition {memneeded > SIZE_MAX}. + if (nfft >= (SIZE_MAX - 2*sizeof(struct kiss_fft_state))/sizeof(kiss_fft_cpx)) + return NULL; + + size_t memneeded = KISS_FFT_ALIGN_SIZE_UP(sizeof(struct kiss_fft_state) + + sizeof(kiss_fft_cpx)*(nfft-1)); /* twiddle factors*/ + + if ( lenmem==NULL ) { + st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded ); + }else{ + if (mem != NULL && *lenmem >= memneeded) + st = (kiss_fft_cfg)mem; + *lenmem = memneeded; + } + if (st) { + int i; + st->nfft=nfft; + st->inverse = inverse_fft; + + for (i=0;iinverse) + phase *= -1; + kf_cexp(st->twiddles+i, phase ); + } + + kf_factor(nfft,st->factors); + } + return st; +} + + +void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride) +{ + if (fin == fout) { + //NOTE: this is not really an in-place FFT algorithm. + //It just performs an out-of-place FFT into a temp buffer + if (fout == NULL){ + KISS_FFT_ERROR("fout buffer NULL."); + return; + } + + kiss_fft_cpx * tmpbuf = (kiss_fft_cpx*)KISS_FFT_TMP_ALLOC( sizeof(kiss_fft_cpx)*st->nfft); + if (tmpbuf == NULL){ + KISS_FFT_ERROR("Memory allocation error."); + return; + } + + + + kf_work(tmpbuf,fin,1,in_stride, st->factors,st); + memcpy(fout,tmpbuf,sizeof(kiss_fft_cpx)*st->nfft); + KISS_FFT_TMP_FREE(tmpbuf); + }else{ + kf_work( fout, fin, 1,in_stride, st->factors,st ); + } +} + +void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) +{ + kiss_fft_stride(cfg,fin,fout,1); +} + + +void kiss_fft_cleanup(void) +{ + // nothing needed any more +} + +int kiss_fft_next_fast_size(int n) +{ + while(1) { + int m=n; + while ( (m%2) == 0 ) m/=2; + while ( (m%3) == 0 ) m/=3; + while ( (m%5) == 0 ) m/=5; + if (m<=1) + break; /* n is completely factorable by twos, threes, and fives */ + n++; + } + return n; +} diff --git a/external/gpr3dviewer/third_party/kissfft/kiss_fft.h b/external/gpr3dviewer/third_party/kissfft/kiss_fft.h new file mode 100644 index 0000000..51b6306 --- /dev/null +++ b/external/gpr3dviewer/third_party/kissfft/kiss_fft.h @@ -0,0 +1,184 @@ +/* + * Copyright (c) 2003-2010, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#ifndef KISS_FFT_H +#define KISS_FFT_H + +#include +#include +#include +#include + +// Define KISS_FFT_SHARED macro to properly export symbols +#ifdef KISS_FFT_SHARED +# ifdef _WIN32 +# ifdef KISS_FFT_BUILD +# define KISS_FFT_API __declspec(dllexport) +# else +# define KISS_FFT_API __declspec(dllimport) +# endif +# else +# define KISS_FFT_API __attribute__ ((visibility ("default"))) +# endif +#else +# define KISS_FFT_API +#endif + +#ifdef __cplusplus +extern "C" { +#endif + +/* + ATTENTION! + If you would like a : + -- a utility that will handle the caching of fft objects + -- real-only (no imaginary time component ) FFT + -- a multi-dimensional FFT + -- a command-line utility to perform ffts + -- a command-line utility to perform fast-convolution filtering + + Then see kfc.h kiss_fftr.h kiss_fftnd.h fftutil.c kiss_fastfir.c + in the tools/ directory. +*/ + +/* User may override KISS_FFT_MALLOC and/or KISS_FFT_FREE. */ +#ifdef USE_SIMD +#ifdef HAVE_LASX +# include +# define kiss_fft_scalar __m256 +# ifndef KISS_FFT_MALLOC +# define KISS_FFT_MALLOC(nbytes) aligned_alloc(32, KISS_FFT_ALIGN_SIZE_UP(nbytes)) +# define KISS_FFT_ALIGN_CHECK(ptr) +# define KISS_FFT_ALIGN_SIZE_UP(size) ((size + 31UL) & ~0x1FUL) +# endif +# ifndef KISS_FFT_FREE +# define KISS_FFT_FREE free +# endif +#elif defined(HAVE_LSX) +# include +# define kiss_fft_scalar __m128 +# ifndef KISS_FFT_MALLOC +# define KISS_FFT_MALLOC(nbytes) aligned_alloc(16, KISS_FFT_ALIGN_SIZE_UP(nbytes)) +# define KISS_FFT_ALIGN_CHECK(ptr) +# define KISS_FFT_ALIGN_SIZE_UP(size) ((size + 15UL) & ~0xFUL) +# endif +# ifndef KISS_FFT_FREE +# define KISS_FFT_FREE free +# endif +#else +# include +# define kiss_fft_scalar __m128 +# ifndef KISS_FFT_MALLOC +# define KISS_FFT_MALLOC(nbytes) _mm_malloc(nbytes,16) +# define KISS_FFT_ALIGN_CHECK(ptr) +# define KISS_FFT_ALIGN_SIZE_UP(size) ((size + 15UL) & ~0xFUL) +# endif +# ifndef KISS_FFT_FREE +# define KISS_FFT_FREE _mm_free +# endif +#endif +#else +# define KISS_FFT_ALIGN_CHECK(ptr) +# define KISS_FFT_ALIGN_SIZE_UP(size) (size) +# ifndef KISS_FFT_MALLOC +# define KISS_FFT_MALLOC malloc +# endif +# ifndef KISS_FFT_FREE +# define KISS_FFT_FREE free +# endif +#endif + + +#ifdef FIXED_POINT +#include +# if (FIXED_POINT == 32) +# define kiss_fft_scalar int32_t +# else +# define kiss_fft_scalar int16_t +# endif +#else +# ifndef kiss_fft_scalar +/* default is float */ +# define kiss_fft_scalar float +# endif +#endif + +typedef struct { + kiss_fft_scalar r; + kiss_fft_scalar i; +}kiss_fft_cpx; + +typedef struct kiss_fft_state* kiss_fft_cfg; + +/* + * kiss_fft_alloc + * + * Initialize a FFT (or IFFT) algorithm's cfg/state buffer. + * + * typical usage: kiss_fft_cfg mycfg=kiss_fft_alloc(1024,0,NULL,NULL); + * + * The return value from fft_alloc is a cfg buffer used internally + * by the fft routine or NULL. + * + * If lenmem is NULL, then kiss_fft_alloc will allocate a cfg buffer using malloc. + * The returned value should be free()d when done to avoid memory leaks. + * + * The state can be placed in a user supplied buffer 'mem': + * If lenmem is not NULL and mem is not NULL and *lenmem is large enough, + * then the function places the cfg in mem and the size used in *lenmem + * and returns mem. + * + * If lenmem is not NULL and ( mem is NULL or *lenmem is not large enough), + * then the function returns NULL and places the minimum cfg + * buffer size in *lenmem. + * */ + +kiss_fft_cfg KISS_FFT_API kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem); + +/* + * kiss_fft(cfg,in_out_buf) + * + * Perform an FFT on a complex input buffer. + * for a forward FFT, + * fin should be f[0] , f[1] , ... ,f[nfft-1] + * fout will be F[0] , F[1] , ... ,F[nfft-1] + * Note that each element is complex and can be accessed like + f[k].r and f[k].i + * */ +void KISS_FFT_API kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout); + +/* + A more generic version of the above function. It reads its input from every Nth sample. + * */ +void KISS_FFT_API kiss_fft_stride(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int fin_stride); + +/* If kiss_fft_alloc allocated a buffer, it is one contiguous + buffer and can be simply free()d when no longer needed*/ +#define kiss_fft_free KISS_FFT_FREE + +/* + Cleans up some memory that gets managed internally. Not necessary to call, but it might clean up + your compiler output to call this before you exit. +*/ +void KISS_FFT_API kiss_fft_cleanup(void); + + +/* + * Returns the smallest integer k, such that k>=n and k has only "fast" factors (2,3,5) + */ +int KISS_FFT_API kiss_fft_next_fast_size(int n); + +/* for real ffts, we need an even size */ +#define kiss_fftr_next_fast_size_real(n) \ + (kiss_fft_next_fast_size( ((n)+1)>>1)<<1) + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/external/gpr3dviewer/third_party/kissfft/kiss_fft_log.h b/external/gpr3dviewer/third_party/kissfft/kiss_fft_log.h new file mode 100644 index 0000000..5012474 --- /dev/null +++ b/external/gpr3dviewer/third_party/kissfft/kiss_fft_log.h @@ -0,0 +1,36 @@ +/* + * Copyright (c) 2003-2010, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#ifndef kiss_fft_log_h +#define kiss_fft_log_h + +#define ERROR 1 +#define WARNING 2 +#define INFO 3 +#define DEBUG 4 + +#define STRINGIFY(x) #x +#define TOSTRING(x) STRINGIFY(x) + +#if defined(NDEBUG) +# define KISS_FFT_LOG_MSG(severity, ...) ((void)0) +#else +# define KISS_FFT_LOG_MSG(severity, ...) \ + fprintf(stderr, "[" #severity "] " __FILE__ ":" TOSTRING(__LINE__) " "); \ + fprintf(stderr, __VA_ARGS__); \ + fprintf(stderr, "\n") +#endif + +#define KISS_FFT_ERROR(...) KISS_FFT_LOG_MSG(ERROR, __VA_ARGS__) +#define KISS_FFT_WARNING(...) KISS_FFT_LOG_MSG(WARNING, __VA_ARGS__) +#define KISS_FFT_INFO(...) KISS_FFT_LOG_MSG(INFO, __VA_ARGS__) +#define KISS_FFT_DEBUG(...) KISS_FFT_LOG_MSG(DEBUG, __VA_ARGS__) + + + +#endif /* kiss_fft_log_h */ diff --git a/external/gpr3dviewer/third_party/kissfft/kiss_fftr.c b/external/gpr3dviewer/third_party/kissfft/kiss_fftr.c new file mode 100644 index 0000000..08b6d7b --- /dev/null +++ b/external/gpr3dviewer/third_party/kissfft/kiss_fftr.c @@ -0,0 +1,171 @@ +/* + * Copyright (c) 2003-2004, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#include "kiss_fftr.h" +#include "_kiss_fft_guts.h" + +struct kiss_fftr_state{ + kiss_fft_cfg substate; + kiss_fft_cpx * tmpbuf; + kiss_fft_cpx * super_twiddles; +#ifdef USE_SIMD + void * pad; +#endif +}; + +kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem) +{ + KISS_FFT_ALIGN_CHECK(mem) + + int i; + kiss_fftr_cfg st = NULL; + size_t subsize = 0, memneeded; + + if (nfft & 1) { + KISS_FFT_ERROR("Real FFT optimization must be even."); + return NULL; + } + nfft >>= 1; + + kiss_fft_alloc (nfft, inverse_fft, NULL, &subsize); + memneeded = sizeof(struct kiss_fftr_state) + subsize + sizeof(kiss_fft_cpx) * ( nfft * 3 / 2); + + if (lenmem == NULL) { + st = (kiss_fftr_cfg) KISS_FFT_MALLOC (memneeded); + } else { + if (*lenmem >= memneeded) + st = (kiss_fftr_cfg) mem; + *lenmem = memneeded; + } + if (!st) + return NULL; + + st->substate = (kiss_fft_cfg) (st + 1); /*just beyond kiss_fftr_state struct */ + st->tmpbuf = (kiss_fft_cpx *) (((char *) st->substate) + subsize); + st->super_twiddles = st->tmpbuf + nfft; + kiss_fft_alloc(nfft, inverse_fft, st->substate, &subsize); + + for (i = 0; i < nfft/2; ++i) { + double phase = + -3.14159265358979323846264338327 * ((double) (i+1) / nfft + .5); + if (inverse_fft) + phase *= -1; + kf_cexp (st->super_twiddles+i,phase); + } + return st; +} + +void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata) +{ + /* input buffer timedata is stored row-wise */ + int k,ncfft; + kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc; + + if ( st->substate->inverse) { + KISS_FFT_ERROR("kiss fft usage error: improper alloc"); + return;/* The caller did not call the correct function */ + } + + ncfft = st->substate->nfft; + + /*perform the parallel fft of two real signals packed in real,imag*/ + kiss_fft( st->substate , (const kiss_fft_cpx*)timedata, st->tmpbuf ); + /* The real part of the DC element of the frequency spectrum in st->tmpbuf + * contains the sum of the even-numbered elements of the input time sequence + * The imag part is the sum of the odd-numbered elements + * + * The sum of tdc.r and tdc.i is the sum of the input time sequence. + * yielding DC of input time sequence + * The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1... + * yielding Nyquist bin of input time sequence + */ + + tdc.r = st->tmpbuf[0].r; + tdc.i = st->tmpbuf[0].i; + C_FIXDIV(tdc,2); + CHECK_OVERFLOW_OP(tdc.r ,+, tdc.i); + CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i); + freqdata[0].r = tdc.r + tdc.i; + freqdata[ncfft].r = tdc.r - tdc.i; +#ifdef USE_SIMD +#ifdef HAVE_LASX + freqdata[0].i = (__m256)(__lasx_xvreplgr2vr_w(0)); + freqdata[ncfft].i = freqdata[0].i; +#elif defined(HAVE_LSX) + freqdata[0].i = (__m128)(__lsx_vreplgr2vr_w(0)); + freqdata[ncfft].i = freqdata[0].i; +#else + freqdata[ncfft].i = freqdata[0].i = _mm_set1_ps(0); +#endif +#else + freqdata[ncfft].i = freqdata[0].i = 0; +#endif + + for ( k=1;k <= ncfft/2 ; ++k ) { + fpk = st->tmpbuf[k]; + fpnk.r = st->tmpbuf[ncfft-k].r; + fpnk.i = - st->tmpbuf[ncfft-k].i; + C_FIXDIV(fpk,2); + C_FIXDIV(fpnk,2); + + C_ADD( f1k, fpk , fpnk ); + C_SUB( f2k, fpk , fpnk ); + C_MUL( tw , f2k , st->super_twiddles[k-1]); + + freqdata[k].r = HALF_OF(f1k.r + tw.r); + freqdata[k].i = HALF_OF(f1k.i + tw.i); + freqdata[ncfft-k].r = HALF_OF(f1k.r - tw.r); + freqdata[ncfft-k].i = HALF_OF(tw.i - f1k.i); + } +} + +void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata) +{ + /* input buffer timedata is stored row-wise */ + int k, ncfft; + + if (st->substate->inverse == 0) { + KISS_FFT_ERROR("kiss fft usage error: improper alloc"); + return;/* The caller did not call the correct function */ + } + + ncfft = st->substate->nfft; + + st->tmpbuf[0].r = freqdata[0].r + freqdata[ncfft].r; + st->tmpbuf[0].i = freqdata[0].r - freqdata[ncfft].r; + C_FIXDIV(st->tmpbuf[0],2); + + for (k = 1; k <= ncfft / 2; ++k) { + kiss_fft_cpx fk, fnkc, fek, fok, tmp; + fk = freqdata[k]; + fnkc.r = freqdata[ncfft - k].r; + fnkc.i = -freqdata[ncfft - k].i; + C_FIXDIV( fk , 2 ); + C_FIXDIV( fnkc , 2 ); + + C_ADD (fek, fk, fnkc); + C_SUB (tmp, fk, fnkc); + C_MUL (fok, tmp, st->super_twiddles[k-1]); + C_ADD (st->tmpbuf[k], fek, fok); + C_SUB (st->tmpbuf[ncfft - k], fek, fok); +#ifdef USE_SIMD +#ifdef HAVE_LASX + __m256 neg_one = (__m256)__lasx_xvreplgr2vr_w(0xBF800000); // -1.0f + st->tmpbuf[ncfft - k].i = __lasx_xvfmul_s(st->tmpbuf[ncfft - k].i, neg_one); +#elif defined(HAVE_LSX) + __m128 neg_one = (__m128)__lsx_vreplgr2vr_w(0xBF800000); // -1.0f + st->tmpbuf[ncfft - k].i = __lsx_vfmul_s(st->tmpbuf[ncfft - k].i, neg_one); +#else + st->tmpbuf[ncfft - k].i *= _mm_set1_ps(-1.0); +#endif +#else + st->tmpbuf[ncfft - k].i *= -1; +#endif + } + kiss_fft (st->substate, st->tmpbuf, (kiss_fft_cpx *) timedata); +} diff --git a/external/gpr3dviewer/third_party/kissfft/kiss_fftr.h b/external/gpr3dviewer/third_party/kissfft/kiss_fftr.h new file mode 100644 index 0000000..7fd73d2 --- /dev/null +++ b/external/gpr3dviewer/third_party/kissfft/kiss_fftr.h @@ -0,0 +1,54 @@ +/* + * Copyright (c) 2003-2004, Mark Borgerding. All rights reserved. + * This file is part of KISS FFT - https://github.com/mborgerding/kissfft + * + * SPDX-License-Identifier: BSD-3-Clause + * See COPYING file for more information. + */ + +#ifndef KISS_FTR_H +#define KISS_FTR_H + +#include "kiss_fft.h" +#ifdef __cplusplus +extern "C" { +#endif + + +/* + + Real optimized version can save about 45% cpu time vs. complex fft of a real seq. + + + + */ + +typedef struct kiss_fftr_state *kiss_fftr_cfg; + + +kiss_fftr_cfg KISS_FFT_API kiss_fftr_alloc(int nfft,int inverse_fft,void * mem, size_t * lenmem); +/* + nfft must be even + + If you don't care to allocate space, use mem = lenmem = NULL +*/ + + +void KISS_FFT_API kiss_fftr(kiss_fftr_cfg cfg,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata); +/* + input timedata has nfft scalar points + output freqdata has nfft/2+1 complex points +*/ + +void KISS_FFT_API kiss_fftri(kiss_fftr_cfg cfg,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata); +/* + input freqdata has nfft/2+1 complex points + output timedata has nfft scalar points +*/ + +#define kiss_fftr_free KISS_FFT_FREE + +#ifdef __cplusplus +} +#endif +#endif diff --git a/tools/gpr3dv_smoke/CMakeLists.txt b/tools/gpr3dv_smoke/CMakeLists.txt new file mode 100644 index 0000000..05f0fd4 --- /dev/null +++ b/tools/gpr3dv_smoke/CMakeLists.txt @@ -0,0 +1,22 @@ +# gpr3dv-smoke —— 冒烟 CLI，验证 vendored 3DGPRViewer 数据生成链编过且跑通。 +# 用法: gpr3dv-smoke +# 例: gpr3dv-smoke "D:/Downloads/明星路" "明星路_001" +# 走原版 API：IprhParser::loadImpulseMultiChannel → GPRDataModel::buildVolumeData +# → RadarProcessor::runPipeline(默认流水线)。打印维度/处理前后统计/通道数。 + +add_executable(gpr3dv_smoke main.cpp) + +set_target_properties(gpr3dv_smoke PROPERTIES + OUTPUT_NAME "gpr3dv-smoke" + AUTOMOC OFF AUTOUIC OFF AUTORCC OFF) + +target_link_libraries(gpr3dv_smoke PRIVATE geopro_gpr3dv Qt6::Core Qt6::Gui) +target_compile_features(gpr3dv_smoke PRIVATE cxx_std_17) + +if(WIN32) + # 运行时 DLL（Qt6Core 等）拷到 exe 旁，使其可直接运行（无需手设 PATH）。 + add_custom_command(TARGET gpr3dv_smoke POST_BUILD + COMMAND ${CMAKE_COMMAND} -E copy_if_different + $ $ + COMMAND_EXPAND_LISTS) +endif() diff --git a/tools/gpr3dv_smoke/main.cpp b/tools/gpr3dv_smoke/main.cpp new file mode 100644 index 0000000..403159a --- /dev/null +++ b/tools/gpr3dv_smoke/main.cpp @@ -0,0 +1,115 @@ +// gpr3dv-smoke —— vendored 3DGPRViewer 数据生成链冒烟入口。 +// +// 验证目标（任务 P1 验收）： +// ① 读一条多通道测线 → GPRDataModel（traces） +// ② GPRDataModel::buildVolumeData → 立方体 volumeData[ch][trace][sample] +// ③ RadarProcessor::runPipeline(默认流水线) → 处理后 GPRDataModel +// 打印：立方体维度 / 处理前后统计量（证明处理生效）/ 通道数（读自数据）。 +// +// 用法: +// gpr3dv-smoke +// 例: gpr3dv-smoke "D:/Downloads/明星路" "明星路_001" +// +// 算法零改动：本文件只调用 vendored 原版静态 API，不触碰任何算法体。 + +#include +#include +#include + +#include + +#include "GPRDataModel.h" +#include "IprhParser.h" +#include "RadarProcessor.h" + +// 全体 traces 的平均绝对幅值——朴素的"处理是否生效"标量证据。 +// RadarProcessor 流水线在 model.traces 上做原位变换（非 volumeData），故统计取自 traces。 +static double meanAbsAmplitude(const GPRDataModel &model, qint64 *sampleCountOut = nullptr) +{ + long double sum = 0.0L; + qint64 count = 0; + for (const RadarTrace &tr : model.traces) { + for (short v : tr.amplitudes) { + sum += static_cast(v < 0 ? -v : v); + ++count; + } + } + if (sampleCountOut) { + *sampleCountOut = count; + } + return count > 0 ? static_cast(sum / count) : 0.0; +} + +int main(int argc, char **argv) +{ + if (argc < 3) { + std::fprintf(stderr, + "用法: gpr3dv-smoke \n" + "例: gpr3dv-smoke \"D:/Downloads/明星路\" \"明星路_001\"\n"); + return 2; + } + + const QString lineDir = QString::fromLocal8Bit(argv[1]); + const QString linePrefix = QString::fromLocal8Bit(argv[2]); + + std::printf("[gpr3dv-smoke] lineDir=%s linePrefix=%s\n", + lineDir.toLocal8Bit().constData(), + linePrefix.toLocal8Bit().constData()); + + // ---- ① 读多通道测线（通道数读自数据，不写死）---- + GPRDataModel model; + if (!IprhParser::loadImpulseMultiChannel(lineDir, linePrefix, model)) { + std::fprintf(stderr, "[gpr3dv-smoke] 错误: loadImpulseMultiChannel 失败\n"); + return 1; + } + + const int channels = model.header.numberOfChannels; // 读自 .iprh CHANNELS / _A 文件数 + const int tracesPerChannel = model.getTracesPerChannel(); + const int samplesPerTrace = model.header.samplesPerTrace; + + std::printf("[gpr3dv-smoke] 加载完成: 通道数=%d 每通道道数=%d 每道样本数=%d 总道数=%lld\n", + channels, tracesPerChannel, samplesPerTrace, + static_cast(model.traces.size())); + + // ---- ② 建三维立方体 volumeData[ch][trace][sample] ---- + model.buildVolumeData(); + const int volCh = model.volumeData.size(); + const int volTr = (volCh > 0) ? model.volumeData[0].size() : 0; + const int volSm = (volCh > 0 && volTr > 0) ? model.volumeData[0][0].size() : 0; + std::printf("[gpr3dv-smoke] 立方体维度 [通道 x 道 x 样本] = %d x %d x %d\n", + volCh, volTr, volSm); + + short gMin = 0, gMax = 0; + model.getGlobalMinMax(gMin, gMax); + std::printf("[gpr3dv-smoke] 立方体全局幅值范围: min=%d max=%d\n", + static_cast(gMin), static_cast(gMax)); + + // ---- 处理前统计 ---- + qint64 nBefore = 0; + const double meanBefore = meanAbsAmplitude(model, &nBefore); + std::printf("[gpr3dv-smoke] 处理前 平均绝对幅值=%.4f (样本数=%lld)\n", + meanBefore, static_cast(nBefore)); + + // ---- ③ RadarProcessor 默认流水线 ---- + RadarProcessor::ProcPipeline pipeline; + pipeline.setDefaultFlow(); + std::printf("[gpr3dv-smoke] 运行默认流水线（步骤数=%lld）...\n", + static_cast(pipeline.steps.size())); + + const GPRDataModel processed = RadarProcessor::runPipeline(model, pipeline); + + qint64 nAfter = 0; + const double meanAfter = meanAbsAmplitude(processed, &nAfter); + std::printf("[gpr3dv-smoke] 处理后 平均绝对幅值=%.4f (样本数=%lld 道数=%lld 样本/道=%d)\n", + meanAfter, static_cast(nAfter), + static_cast(processed.traces.size()), + processed.header.samplesPerTrace); + + const double delta = meanAfter - meanBefore; + std::printf("[gpr3dv-smoke] 处理前后差值=%.4f —— %s\n", + delta, + (qAbs(delta) > 1e-9 ? "处理已生效(统计量变化)" : "警告:统计量未变化")); + + std::printf("[gpr3dv-smoke] OK\n"); + return 0; +}