sunwen/UR
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add transmossion detection functions.

Browse Source
This commit is contained in:
kradchen
2023-05-23 09:39:21 +08:00
parent b84eb97bbe
commit 6643752cad
4 changed files with 347 additions and 68 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
CMakeLists.txt
View File

@@ -6,6 +6,7 @@ set(CMAKE_INCLUDE_CURRENT_DIR ON)
find_package(Aurora REQUIRED)
find_package(Parser REQUIRED)
find_package(URDepends REQUIRED)
file(GLOB_RECURSE cpp_files ./src/*.cpp)
file(GLOB_RECURSE cxx_files ./src/*.cxx)
file(GLOB_RECURSE header_files ./src/*.h)
@@ -13,9 +14,11 @@ add_executable(UR ${cpp_files} ${cxx_files} ${header_files} ${Aurora_Source})
target_compile_options(UR PUBLIC ${Aurora_Complie_Options} "-march=native")
target_include_directories(UR PUBLIC ${Aurora_INCLUDE_DIRS})
target_include_directories(UR PUBLIC ${Parser_INCLUDE_DIRS})
target_include_directories(UR PUBLIC ${URDepends_INCLUDES_DIRS})
target_link_libraries(UR PUBLIC ${Aurora_Libraries})
target_link_libraries(UR PUBLIC matio)
target_link_libraries(UR PUBLIC ${Parser_Libraries})
target_link_libraries(UR PUBLIC URDepends::TransDetection)
find_package(GTest REQUIRED)
INCLUDE_DIRECTORIES(${GTEST_INCLUDE_DIRS})
file(GLOB_RECURSE test_cpp ./test/*.cpp)
@@ -26,8 +29,10 @@ target_include_directories(UR_Test PUBLIC ./test/ ./src/)
target_compile_options(UR_Test PUBLIC ${Aurora_Complie_Options} "-march=native")
target_include_directories(UR_Test PUBLIC ${Aurora_INCLUDE_DIRS})
target_include_directories(UR_Test PUBLIC ${Parser_INCLUDE_DIRS})
target_include_directories(UR_Test PUBLIC ${URDepends_INCLUDES_DIRS})
target_link_libraries(UR_Test PUBLIC ${Aurora_Libraries})
target_link_libraries(UR_Test PUBLIC matio)
target_link_libraries(UR_Test PUBLIC ${Parser_Libraries})
target_link_libraries(UR_Test PUBLIC URDepends::TransDetection)
target_link_libraries(UR_Test PUBLIC ${GTEST_BOTH_LIBRARIES})
gtest_discover_tests(UR_Test )

246
src/transmissionReconstruction/detection/detection.cpp
View File

@@ -1,10 +1,16 @@
#include "detection.h"
#include <algorithm>
#include <cstddef>
#include <sys/types.h>
#include "Function.h"
#include "Function1D.h"
#include "Function2D.h"
#include "Function3D.h"
#include "Matrix.h"
#include <sys/types.h>
#include "calculateBankDetectAndHilbertTransformation.hpp"
using namespace Aurora;
namespace Recon {
@@ -141,10 +147,10 @@ namespace Recon {
auto endPosRef = zeros(Ascan.getDimSize(1), 1);
for (size_t i = 0; i < Ascan.getDimSize(1); i++)
{
startPos(i) = floor(max(tof(i).toMatrix()*sampleRate+offsetElectronicSamples,1));
endPos(i) = ceil(min(sizeAscan(1).toMatrix(), tof(i).toMatrix()*sampleRate+offsetElectronicSamples+detectionWindowATT));
startPosRef(i) = floor(max( tof2(i).toMatrix(),1));
endPosRef(i) = ceil(min(sizeAscan(1).toMatrix(), tof2(i).toMatrix()+detectionWindowATT));
startPos[i] = std::floor(std::max(tof[i]*sampleRate+offsetElectronicSamples,1.0));
endPos[i] = std::ceil(std::min(sizeAscan[1], tof[i]*sampleRate+offsetElectronicSamples+detectionWindowATT));
startPosRef[i] = std::floor(std::max( tof2[i],1.0));
endPosRef[i] = std::ceil(std::min(sizeAscan[1], tof2[i]+detectionWindowATT));
}
return calculateAttenuation(envelopeOfAScan,startPos,endPos,envelopeOfReferenceAScan,startPosRef,endPosRef);
}
@@ -158,8 +164,8 @@ namespace Recon {
}
return result;
}
//TODO:need test
Aurora::Matrix detectTofVectorized(
DetectResult detectTofVectorized(
const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
const Aurora::Matrix &distBlock, const Aurora::Matrix &distBlockRef,
const Aurora::Matrix &sosWaterBlock,
@@ -175,10 +181,15 @@ namespace Recon {
timeResult1.AscanBlockProcessed = AscanBlock;
TimeWindowResult timeResult2;
timeResult2.AscanBlockProcessed = AscanRefBlock;
if (useTimeWindowing == 1)
{
timeResult1 = applyTimeWindowing(AscanBlock, sampleRate, distBlock, sosWaterBlock, expectedSOSWater, offsetElectronicSamples, detectionWindowSOS, minSpeedOfSound, maxSpeedOfSound, gaussWindow);
timeResult2 = applyTimeWindowing(AscanRefBlock, sampleRate, distBlockRef, sosWaterRefBlock, expectedSOSWater, offsetElectronicSamples, detectionWindowSOS, minSpeedOfSound, maxSpeedOfSound, gaussWindow);
if (useTimeWindowing == 1) {
timeResult1 = applyTimeWindowing(
AscanBlock, sampleRate, distBlock, sosWaterBlock,
expectedSOSWater, offsetElectronicSamples, detectionWindowSOS,
minSpeedOfSound, maxSpeedOfSound, gaussWindow);
timeResult2 = applyTimeWindowing(
AscanRefBlock, sampleRate, distBlockRef, sosWaterRefBlock,
expectedSOSWater, offsetElectronicSamples, detectionWindowSOS,
minSpeedOfSound, maxSpeedOfSound, gaussWindow);
diffStartSearch = timeResult1.startSearch - timeResult1.startSearch;
}
@@ -186,40 +197,213 @@ namespace Recon {
auto _AscanRefBlock = timeResult2.AscanBlockProcessed;
auto m = std::max(size(_AscanBlock, 1), size(_AscanRefBlock, 1));
auto maxlag = std::max(size(_AscanBlock, 1), size(_AscanRefBlock, 1)) - 1;
auto mxl = std::min(maxlag, m-1);
auto ceilLog2 = nextpow2(2*m-1);
auto m2 = pow(2,ceilLog2);
auto maxlag =
std::max(size(_AscanBlock, 1), size(_AscanRefBlock, 1)) - 1;
auto mxl = std::min(maxlag, m - 1);
auto ceilLog2 = nextpow2(2 * m - 1);
auto m2 = pow(2, ceilLog2);
auto x = fft(_AscanBlock, m2);
auto y = fft(_AscanRefBlock, m2);
auto c_1_1 = x*conj(y);
auto c_1_1 = x * conj(y);
auto c_1_2 = ifft(c_1_1);
auto c1 = real(c_1_2);
auto c = zeros(mxl+mxl+1,c1.getDimSize(1));
auto c = zeros(mxl + mxl + 1, c1.getDimSize(1));
#pragma omp parallel for
for (size_t i = 0; i < mxl; i++)
{
c(i,$) = c1(m2-mxl+i, $);
c(i+mxl,$) = c1(i, $);
for (size_t i = 0; i < mxl; i++) {
c(i, $) = c1(m2 - mxl + i, $);
c(i + mxl, $) = c1(i, $);
}
c(mxl+mxl,$) = c1(mxl, $);
c(mxl + mxl, $) = c1(mxl, $);
auto shiftInSamples =zeros(1,c1.getDimSize(1));
auto shiftInSamples = zeros(1, c1.getDimSize(1));
#pragma omp parallel for
for (size_t i = 0; i < c1.getDimSize(1); i++)
{
long rowID=0,colID=0;
max(c($,i).toMatrix(),All,rowID,colID);
shiftInSamples[i]=-maxlag+colID+rowID;
for (size_t i = 0; i < c1.getDimSize(1); i++) {
long rowID = 0, colID = 0;
max(c($, i).toMatrix(), All, rowID, colID);
shiftInSamples[i] = -maxlag + colID + rowID;
}
if (useTimeWindowing)
{
if (useTimeWindowing) {
shiftInSamples = shiftInSamples - diffStartSearch;
}
auto tof = shiftInSamples / sampleRate + distBlock / sosWaterBlock;
auto sosValue = distBlock / tof;
return tof;
DetectResult result;
result.tof = tof;
result.sosValue = sosValue;
return result;
}
//TODO:未测
DetectResult detectTofAndAtt(
const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
const Aurora::Matrix &distBlock, const Aurora::Matrix &distRefBlock,
const Aurora::Matrix &sosWaterBlock,
const Aurora::Matrix &sosWaterRefBlock,
int resampleFactor,int nthreads, double expectedSOSWater,
int useTimeWindowing, int aScanReconstructionFrequency,int detectionWindowATT,
double offsetElectronic, int detectionWindowSOS, double minSpeedOfSound,
double maxSpeedOfSound, bool gaussWindow)
{
auto result = detectTofVectorized(
AscanBlock, AscanRefBlock, distBlock, distRefBlock, sosWaterBlock,
sosWaterRefBlock, expectedSOSWater, useTimeWindowing,
aScanReconstructionFrequency, offsetElectronic, detectionWindowSOS,
minSpeedOfSound, maxSpeedOfSound, gaussWindow);
auto tofRel = result.tof - distBlock / sosWaterBlock;
result.att = detectAttVectorized(
AscanBlock, AscanRefBlock, distRefBlock, sosWaterRefBlock,
result.tof, aScanReconstructionFrequency, offsetElectronic,
detectionWindowATT);
result.tof = tofRel;
return result;
}
//TODO:未测
DetectResult detectTofAndAttMex(
const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
const Aurora::Matrix &distBlock, const Aurora::Matrix &distRefBlock,
const Aurora::Matrix &sosWaterBlock,
const Aurora::Matrix &sosWaterRefBlock,
int resampleFactor,int nthreads, double expectedSOSWater,
int useTimeWindowing, int aScanReconstructionFrequency,int detectionWindowATT,
double offsetElectronic, int detectionWindowSOS, double minSpeedOfSound,
double maxSpeedOfSound, bool gaussWindow)
{
auto sizeAscan = size(AscanBlock);
auto sampleRate = aScanReconstructionFrequency;
double offsetElectronicSamples = offsetElectronic * sampleRate;
Matrix diffStartSearch;
TimeWindowResult timeResult1;
timeResult1.AscanBlockProcessed = AscanBlock;
TimeWindowResult timeResult2;
timeResult2.AscanBlockProcessed = AscanRefBlock;
if (useTimeWindowing == 1) {
timeResult1 = applyTimeWindowing(
AscanBlock, sampleRate, distBlock, sosWaterBlock,
expectedSOSWater, offsetElectronicSamples, detectionWindowSOS,
minSpeedOfSound, maxSpeedOfSound, gaussWindow);
timeResult2 = applyTimeWindowing(
AscanRefBlock, sampleRate, distRefBlock, sosWaterRefBlock,
expectedSOSWater, offsetElectronicSamples, detectionWindowSOS,
minSpeedOfSound, maxSpeedOfSound, gaussWindow);
diffStartSearch = timeResult1.startSearch - timeResult1.startSearch;
}
auto _AscanBlock = timeResult1.AscanBlockProcessed;
auto _AscanRefBlock = timeResult2.AscanBlockProcessed;
int M =std::min(AscanBlock.getDimSize(0),AscanRefBlock.getDimSize(0));
double * resDetectD = nullptr;
double * resEnvelopeD = nullptr;
double * resEnvelopeRefD = nullptr;
size_t N = _AscanBlock.getDimSize(1);
size_t totalSize = N*M;
{
auto resDetect = new float[N]{0};
auto resEnvelope = new float[totalSize]{0};
auto resEnvelopeRef = new float[totalSize]{0};
auto aScans_r = new float[totalSize]{0};
auto aScansRef_r = new float[totalSize]{0};
{
Matrix _AscanBlock_trim = _AscanBlock.getDimSize(0)!=M?_AscanBlock.block(0, 0, M-1):_AscanBlock;
Matrix _AscanRefBlock_trim = _AscanRefBlock.getDimSize(0)!=M?_AscanRefBlock.block(0, 0, M-1):_AscanRefBlock;
std::copy(_AscanBlock_trim.getData(), _AscanBlock_trim.getData() + totalSize, aScans_r);
std::copy(_AscanRefBlock_trim.getData(), _AscanRefBlock_trim.getData() + totalSize,
aScans_r);
}
calculateBankDetectAndHilbertTransformation(
aScans_r, aScansRef_r, M, N, resampleFactor, nthreads,
resDetect, resEnvelope, resEnvelopeRef);
delete [] aScans_r;
delete [] aScansRef_r;
resDetectD = Aurora::malloc(N);
std::copy(resDetect, resDetect+ N, resDetectD);
delete [] resDetect;
resEnvelopeD = Aurora::malloc(totalSize);
std::copy(resEnvelope, resEnvelope+ totalSize, resEnvelopeD);
delete [] resEnvelope;
resEnvelopeRefD = Aurora::malloc(totalSize);
std::copy(resEnvelopeRef, resEnvelopeRef+ totalSize, resEnvelopeRefD);
delete [] resEnvelopeRef;
}
auto resDetect =Matrix::New(resDetectD,1,N);
auto resEnvelope =Matrix::New(resEnvelopeD,N,M);
auto resEnvelopeRef =Matrix::New(resEnvelopeRefD,N,M);
//floor(size(AscanBlock, 1)*inits.resampleFactor / 2 - 1),
int end_1 = std::floor(_AscanBlock.getDimSize(0)*resampleFactor/2-1);
//-ceil(size(AscanBlock, 1)*inits.resampleFactor / 2)
int begin_2 = -std::ceil(_AscanBlock.getDimSize(0)*resampleFactor/2);
int *lags2 = new int[_AscanBlock.getDimSize(0)];
for (size_t i = 0; i <= end_1; i++)
{
lags2[i] = i;
lags2[i+end_1+1] = begin_2+i;
}
auto resDetectLags = zeros(_AscanBlock.getDimSize(1),1);
for (size_t i = 0; i < _AscanBlock.getDimSize(1); i++)
{
resDetectLags[i] = lags2[(int)resDetect[i]+1];
}
delete [] lags2;
resDetectLags =resDetectLags/resampleFactor;
if (useTimeWindowing == 1) {
resDetectLags = resDetectLags - diffStartSearch;
}
auto tofRel = (resDetectLags / sampleRate);
auto tofAbs = tofRel + (distBlock / sosWaterBlock);
auto sosValue = distBlock / tofAbs;
auto tof2 = (distRefBlock / sosWaterRefBlock) * sampleRate + offsetElectronicSamples;
auto startPos = zeros(_AscanBlock.getDimSize(1), 1);
auto endPos = zeros(_AscanBlock.getDimSize(1), 1);
auto startPosRef = zeros(_AscanBlock.getDimSize(1), 1);
auto endPosRef = zeros(_AscanBlock.getDimSize(1), 1);
for (size_t i = 0; i < _AscanBlock.getDimSize(1); i++)
{
startPos[i] = std::floor(std::max(tofAbs[i]*sampleRate+offsetElectronicSamples,1.0));
endPos[i] = std::ceil(std::min(sizeAscan[1], tofAbs[i]*sampleRate+offsetElectronicSamples+detectionWindowATT));
startPosRef[i] = std::floor(std::max( tof2[i],1.0));
endPosRef[i] = std::ceil(std::min(sizeAscan[1], tof2[i]+detectionWindowATT));
}
DetectResult result;
result.att = calculateAttenuation(resEnvelope,startPos,endPos,resEnvelopeRef,startPosRef,endPosRef);
result.tof = tofRel;
result.sosValue = sosValue;
return result;
}
DetectResult transmissionDetection(const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
const Aurora::Matrix &distBlock, const Aurora::Matrix &distRefBlock,
const Aurora::Matrix &sosWaterBlock,
const Aurora::Matrix &sosWaterRefBlock,int version, int resampleFactor,int nthreads, double expectedSOSWater,
int useTimeWindowing, int aScanReconstructionFrequency,int detectionWindowATT,
double offsetElectronic, int detectionWindowSOS, double minSpeedOfSound,
double maxSpeedOfSound, bool gaussWindow){
switch (version) {
case 1: {
return detectTofAndAttMex(
AscanBlock, AscanRefBlock, distBlock, distRefBlock,
sosWaterBlock, sosWaterRefBlock, resampleFactor, nthreads,
expectedSOSWater, useTimeWindowing,
aScanReconstructionFrequency, detectionWindowATT,
offsetElectronic, detectionWindowSOS, minSpeedOfSound,
maxSpeedOfSound, gaussWindow);
}
case 2:
default:
return detectTofAndAtt(
AscanBlock, AscanRefBlock, distBlock, distRefBlock,
sosWaterBlock, sosWaterRefBlock, resampleFactor, nthreads,
expectedSOSWater, useTimeWindowing,
aScanReconstructionFrequency, detectionWindowATT,
offsetElectronic, detectionWindowSOS, minSpeedOfSound,
maxSpeedOfSound, gaussWindow);
}
}
}

41
src/transmissionReconstruction/detection/detection.h
View File

@@ -11,9 +11,10 @@ struct TimeWindowResult {
Aurora::Matrix AscanBlockProcessed;
};
struct DetectTofResult {
struct DetectResult {
Aurora::Matrix tof;
Aurora::Matrix sosValue;
Aurora::Matrix att;
};
Aurora::Matrix calculateAttenuation(const Aurora::Matrix &ascans,
const Aurora::Matrix &startPos,
@@ -35,20 +36,52 @@ TimeWindowResult applyTimeWindowing(
double expectedSOSWater, double startOffset, double segmentLenOffset,
double minSpeedOfSound, double maxSpeedOfSound, bool gaussWindow);
Aurora::Matrix detectAttVectorized(const Aurora::Matrix &Ascan, const Aurora::Matrix &AscanRef,
Aurora::Matrix
detectAttVectorized(const Aurora::Matrix &Ascan, const Aurora::Matrix &AscanRef,
const Aurora::Matrix &distRef,
const Aurora::Matrix &sosWaterRef,
const Aurora::Matrix &tof, int aScanReconstructionFrequency,
double offsetElectronic, int detectionWindowATT);
Aurora::Matrix detectTofVectorized(
DetectResult
detectTofVectorized(
const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
const Aurora::Matrix &distBlock, const Aurora::Matrix &distBlockRef,
const Aurora::Matrix &sosWaterBlock,
const Aurora::Matrix &sosWaterRefBlock, double expectedSOSWater,
int useTimeWindowing, int aScanReconstructionFrequency,
double offsetElectronic, int detectionWindowSOS, double minSpeedOfSound,
double maxSpeedOfSound, bool gaussWindow) ;
double maxSpeedOfSound, bool gaussWindow);
DetectResult
detectTofAndAtt(
const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
const Aurora::Matrix &distBlock, const Aurora::Matrix &distRefBlock,
const Aurora::Matrix &sosWaterBlock, const Aurora::Matrix &sosWaterRefBlock,
int resampleFactor, int nthreads, double expectedSOSWater,
int useTimeWindowing, int aScanReconstructionFrequency,
int detectionWindowATT, double offsetElectronic, int detectionWindowSOS,
double minSpeedOfSound, double maxSpeedOfSound, bool gaussWindow);
DetectResult
detectTofAndAttMex(
const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
const Aurora::Matrix &distBlock, const Aurora::Matrix &distBlockRef,
const Aurora::Matrix &sosWaterBlock, const Aurora::Matrix &sosWaterRefBlock,
int resampleFactor, int nthreads, double expectedSOSWater,
int useTimeWindowing, int aScanReconstructionFrequency,
int detectionWindowATT, double offsetElectronic, int detectionWindowSOS,
double minSpeedOfSound, double maxSpeedOfSound, bool gaussWindow);
DetectResult
transmissionDetection(
const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
const Aurora::Matrix &distBlock, const Aurora::Matrix &distRefBlock,
const Aurora::Matrix &sosWaterBlock, const Aurora::Matrix &sosWaterRefBlock,
int version, int resampleFactor, int nthreads, double expectedSOSWater,
int useTimeWindowing, int aScanReconstructionFrequency,
int detectionWindowATT, double offsetElectronic, int detectionWindowSOS,
double minSpeedOfSound, double maxSpeedOfSound, bool gaussWindow);
} // namespace Recon

57
test/Detection_Test.cpp
View File

@@ -1,3 +1,4 @@
#include <cmath>
#include <gtest/gtest.h>
#include <limits>
@@ -36,6 +37,62 @@ protected:
}
};
TEST_F(Detection_Test, detectTofAndAttMex) {
MatlabReader m("/home/krad/TestData/getBlockOfTransmissionData.mat");
auto AscanBlock = m.read("AscanBlock");
auto AscanRefBlock = m.read("AscanRefBlock");
auto distBlock = m.read("dists");
auto distBlockRef = m.read("distRefBlock");
auto sosWaterBlock = m.read("waterTempBlock");
auto sosWaterRefBlock = m.read("waterTempRefBlock");
double expectedSOSWater = 1.511948131508464e+03;
auto result = Recon::detectTofAndAttMex(
AscanBlock, AscanRefBlock, distBlock, distBlockRef, sosWaterBlock,
sosWaterRefBlock, Recon::transParams::resampleFactor,
Recon::transParams::nThreads, expectedSOSWater,
Recon::transParams::useTimeWindowing,
Recon::transParams::aScanReconstructionFrequency,
Recon::transParams::offsetElectronic,Recon::transParams::detectionWindowATT,
Recon::transParams::detectionWindowSOS,
Recon::transParams::minSpeedOfSound,
Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow);
MatlabReader m2("/home/krad/TestData/tofResult.mat");
auto tof = m2.read("tof");
auto sosvalue = m2.read("sosValue");
EXPECT_EQ(tof.getDataSize(), result.tof.getDataSize());
}
TEST_F(Detection_Test, detectAttVectorized) {
MatlabReader m("/home/krad/TestData/getBlockOfTransmissionData.mat");
auto AscanBlock = m.read("AscanBlock");
auto AscanRefBlock = m.read("AscanRefBlock");
auto distBlockRef = m.read("distRefBlock");
auto sosWaterRefBlock = m.read("waterTempRefBlock");
MatlabReader m2("/home/krad/TestData/tofResult.mat");
auto tof = m2.read("tof");
double expectedSOSWater = 1.511948131508464e+03;
auto result = Recon::detectAttVectorized(
AscanBlock, AscanRefBlock, distBlockRef,sosWaterRefBlock,
tof,
Recon::transParams::aScanReconstructionFrequency,
Recon::transParams::offsetElectronic,
Recon::transParams::detectionWindowSOS);
for (size_t i = 0; i < result.getDataSize(); i++)
{
EXPECT_TRUE(std::isnan(result[i]))<<",index:"<<i;
}
}
TEST_F(Detection_Test, calculateStarEndSearchPosition) {
auto distBlock = Aurora::Matrix::fromRawData(new double[3]{0.22, 0.21, 0.11}, 3, 1);