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

Fix detectTofAndAtt and main log time

Browse Source
This commit is contained in:
kradchen
2023-06-05 14:55:53 +08:00
parent 90fd2afe7c
commit 15cfc29e6f
3 changed files with 114 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
src/main.cxx
View File

@@ -9,8 +9,26 @@
#include "common/ceMatchedFilterHandling.h" #include "common/ceMatchedFilterHandling.h"
#include "MatlabReader.h" #include "MatlabReader.h"
#include "startReconstructions.h" #include "startReconstructions.h"
#include "spdlog/spdlog.h"
#include "spdlog/sinks/stdout_color_sinks.h"
#include "spdlog/sinks/basic_file_sink.h"
std::shared_ptr<spdlog::logger> getLogger(const char* title) {
auto console_sink = std::make_shared<spdlog::sinks::stdout_color_sink_mt>();
console_sink->set_level(spdlog::level::info);
console_sink->set_pattern(fmt::format("[%Y-%m-%d %T .%e][{}] [%^%l%$] %v", title));
std::shared_ptr<spdlog::logger> logger(new spdlog::logger(title, {console_sink}));
logger->set_level(spdlog::level::info);
logger->flush_on(spdlog::level::info);
return logger;
}
int main() int main()
{ {
auto defaultLogger = getLogger("Main");
spdlog::set_default_logger(defaultLogger);
SPDLOG_INFO("start");
Recon::startReconstructions(); Recon::startReconstructions();
SPDLOG_INFO("finish");
return 0; return 0;
} }

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

@@ -151,7 +151,7 @@ namespace Recon {
for (size_t i = 0; i < Ascan.getDimSize(1); i++) for (size_t i = 0; i < Ascan.getDimSize(1); i++)
{ {
startPos[i] = std::floor(std::max(tof[i]*sampleRate+offsetElectronicSamples,1.0)); 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)); endPos[i] = std::ceil(std::min(sizeAscan[0], tof[i]*sampleRate+offsetElectronicSamples+detectionWindowATT));
startPosRef[i] = std::floor(std::max( tof2[i],1.0)); startPosRef[i] = std::floor(std::max( tof2[i],1.0));
endPosRef[i] = std::ceil(std::min(sizeAscan[1], tof2[i]+detectionWindowATT)); endPosRef[i] = std::ceil(std::min(sizeAscan[1], tof2[i]+detectionWindowATT));
} }
@@ -194,7 +194,7 @@ namespace Recon {
expectedSOSWater, offsetElectronicSamples, detectionWindowSOS, expectedSOSWater, offsetElectronicSamples, detectionWindowSOS,
minSpeedOfSound, maxSpeedOfSound, gaussWindow); minSpeedOfSound, maxSpeedOfSound, gaussWindow);
diffStartSearch = timeResult1.startSearch - timeResult1.startSearch; diffStartSearch = timeResult1.startSearch - timeResult2.startSearch;
} }
auto _AscanBlock = timeResult1.AscanBlockProcessed; auto _AscanBlock = timeResult1.AscanBlockProcessed;
auto _AscanRefBlock = timeResult2.AscanBlockProcessed; auto _AscanRefBlock = timeResult2.AscanBlockProcessed;
@@ -213,15 +213,17 @@ namespace Recon {
auto c1 = real(c_1_2); 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 c.setBlock(0, 0, mxl-1, c1.block(0, m2-mxl,m2-1));
for (size_t i = 0; i < mxl; i++) { c.setBlock(0, mxl, mxl*2, c1.block(0, 0, mxl));
c(i, $) = c1(m2 - mxl + i, $); // #pragma omp parallel for
c(i + mxl, $) = c1(i, $); // for (size_t i = 0; i < mxl; i++) {
} // c(i, $) = c1(m2 - mxl + i, $);
c(mxl + mxl, $) = c1(mxl, $); // c(i + mxl, $) = c1(i, $);
// }
// c(mxl + mxl, $) = c1(mxl, $);
auto shiftInSamples = zeros(1, c1.getDimSize(1)); auto shiftInSamples = zeros(1, c1.getDimSize(1));
#pragma omp parallel for // #pragma omp parallel for
for (size_t i = 0; i < c1.getDimSize(1); i++) { for (size_t i = 0; i < c1.getDimSize(1); i++) {
long rowID = 0, colID = 0; long rowID = 0, colID = 0;
max(c($, i).toMatrix(), All, rowID, colID); max(c($, i).toMatrix(), All, rowID, colID);
@@ -248,15 +250,69 @@ namespace Recon {
double offsetElectronic, int detectionWindowSOS, double minSpeedOfSound, double offsetElectronic, int detectionWindowSOS, double minSpeedOfSound,
double maxSpeedOfSound, bool gaussWindow) double maxSpeedOfSound, bool gaussWindow)
{ {
auto result = detectTofVectorized( auto sampleRate = aScanReconstructionFrequency;
AscanBlock, AscanRefBlock, distBlock, distRefBlock, sosWaterBlock, double offsetElectronicSamples = offsetElectronic * sampleRate;
sosWaterRefBlock, expectedSOSWater, useTimeWindowing, Matrix diffStartSearch;
aScanReconstructionFrequency, offsetElectronic, detectionWindowSOS, TimeWindowResult timeResult1;
minSpeedOfSound, maxSpeedOfSound, gaussWindow); timeResult1.AscanBlockProcessed = AscanBlock;
auto tofRel = result.tof - distBlock / sosWaterBlock; 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 - timeResult2.startSearch;
}
auto _AscanBlock = timeResult1.AscanBlockProcessed;
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 x = fft(_AscanBlock, m2);
auto y = fft(_AscanRefBlock, m2);
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));
c.setBlock(0, 0, mxl-1, c1.block(0, m2-mxl,m2-1));
c.setBlock(0, mxl, mxl*2, c1.block(0, 0, mxl));
// #pragma omp parallel for
// for (size_t i = 0; i < mxl; i++) {
// c(i, $) = c1(m2 - mxl + i, $);
// c(i + mxl, $) = c1(i, $);
// }
// c(mxl + mxl, $) = c1(mxl, $);
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;
}
if (useTimeWindowing) {
shiftInSamples = shiftInSamples - diffStartSearch;
}
auto tof = shiftInSamples / sampleRate + distBlock / sosWaterBlock;
auto sosValue = distBlock / tof;
DetectResult result;
result.sosValue = sosValue;
auto tofRel = tof - distBlock / sosWaterBlock;
result.att = detectAttVectorized( result.att = detectAttVectorized(
AscanBlock, AscanRefBlock, distRefBlock, sosWaterRefBlock, _AscanBlock, _AscanRefBlock, distRefBlock, sosWaterRefBlock,
result.tof, aScanReconstructionFrequency, offsetElectronic, tof, aScanReconstructionFrequency, offsetElectronic,
detectionWindowATT); detectionWindowATT);
result.tof = tofRel; result.tof = tofRel;
return result; return result;

30
test/Detection_Test.cpp
View File

@@ -63,22 +63,23 @@ TEST_F(Detection_Test, detectTofAndAttMex) {
MatlabReader m3("/home/krad/TestData/sosResult.mat"); MatlabReader m3("/home/krad/TestData/sosResult.mat");
auto sosvalue = m3.read("sosValue"); auto sosvalue = m3.read("sosValue");
auto tof = m3.read("tofRel");
ASSERT_EQ(sosvalue.getDataSize(), result.tof.getDataSize()); ASSERT_EQ(sosvalue.getDataSize(), result.tof.getDataSize());
ASSERT_EQ(sosvalue.getDataSize(), result.sosValue.getDataSize()); ASSERT_EQ(sosvalue.getDataSize(), result.sosValue.getDataSize());
ASSERT_EQ(sosvalue.getDataSize(), result.att.getDataSize()); ASSERT_EQ(sosvalue.getDataSize(), result.att.getDataSize());
#pragma omp parallel for #pragma omp parallel for
for (size_t i = 0; i < result.tof.getDataSize(); i++) for (size_t i = 0; i < result.tof.getDataSize(); i++)
{ {
EXPECT_DOUBLE_AE(0,result.tof[i])<<",index:"<<i; EXPECT_DOUBLE_AE(tof[i],result.tof[i])<<",index:"<<i;
EXPECT_DOUBLE_AE(sosvalue[i],result.sosValue[i])<<",index:"<<i; EXPECT_DOUBLE_AE(sosvalue[i],result.sosValue[i])<<",index:"<<i;
EXPECT_TRUE(result.att[i]==0)<<",index:"<<i; // EXPECT_TRUE(result.att[i]==0)<<",index:"<<i;
} }
} }
TEST_F(Detection_Test, detectTofAndAtt) { TEST_F(Detection_Test, detectTofAndAtt) {
MatlabReader m("/home/krad/TestData/getBlockOfTransmissionData.mat"); MatlabReader m("/home/sun/testData/transmissionDetection.mat");
auto AscanBlock = m.read("AscanBlock"); auto AscanBlock = m.read("AscanBlock");
auto AscanRefBlock = m.read("AscanRefBlock"); auto AscanRefBlock = m.read("AscanRefBlock");
@@ -86,7 +87,8 @@ TEST_F(Detection_Test, detectTofAndAtt) {
auto distBlockRef = m.read("distRefBlock"); auto distBlockRef = m.read("distRefBlock");
auto sosWaterBlock = Recon::temperatureToSoundSpeed(m.read("waterTempBlock"), "marczak"); auto sosWaterBlock = Recon::temperatureToSoundSpeed(m.read("waterTempBlock"), "marczak");
auto sosWaterRefBlock = Recon::temperatureToSoundSpeed(m.read("waterTempRefBlock"), "marczak"); auto sosWaterRefBlock = Recon::temperatureToSoundSpeed(m.read("waterTempRefBlock"), "marczak");
double expectedSOSWater = 1.511948131508464e+03; double expectedSOSWater = 1.512677498767504e+03;
auto result = Recon::detectTofAndAtt( auto result = Recon::detectTofAndAtt(
AscanBlock, AscanRefBlock, distBlock, distBlockRef, sosWaterBlock, AscanBlock, AscanRefBlock, distBlock, distBlockRef, sosWaterBlock,
@@ -94,22 +96,36 @@ TEST_F(Detection_Test, detectTofAndAtt) {
Recon::transParams::nThreads, expectedSOSWater, Recon::transParams::nThreads, expectedSOSWater,
Recon::transParams::useTimeWindowing, Recon::transParams::useTimeWindowing,
Recon::transParams::aScanReconstructionFrequency, Recon::transParams::aScanReconstructionFrequency,
Recon::transParams::offsetElectronic,Recon::transParams::detectionWindowATT, Recon::transParams::detectionWindowATT,Recon::transParams::offsetElectronic,
Recon::transParams::detectionWindowSOS, Recon::transParams::detectionWindowSOS,
Recon::transParams::minSpeedOfSound, Recon::transParams::minSpeedOfSound,
Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow); Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow);
MatlabReader m2("/home/krad/TestData/sosResult.mat"); MatlabReader m2("/home/krad/TestData/sosResult.mat");
auto sosvalue = m2.read("sosValue"); auto sosvalue = m2.read("sosValue");
auto tof = m2.read("tofRel");
auto att = m2.read("att");
// auto result1 = Recon::detectTofAndAttMex(
// AscanBlock, AscanRefBlock, distBlock, distBlockRef, sosWaterBlock,
// sosWaterRefBlock, Recon::transParams::resampleFactor,
// Recon::transParams::nThreads, expectedSOSWater,
// Recon::transParams::useTimeWindowing,
// Recon::transParams::aScanReconstructionFrequency,
// Recon::transParams::detectionWindowATT,Recon::transParams::offsetElectronic,
// Recon::transParams::detectionWindowSOS,
// Recon::transParams::minSpeedOfSound,
// Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow);
ASSERT_EQ(sosvalue.getDataSize(), result.tof.getDataSize()); ASSERT_EQ(sosvalue.getDataSize(), result.tof.getDataSize());
ASSERT_EQ(sosvalue.getDataSize(), result.sosValue.getDataSize()); ASSERT_EQ(sosvalue.getDataSize(), result.sosValue.getDataSize());
ASSERT_EQ(sosvalue.getDataSize(), result.att.getDataSize()); ASSERT_EQ(sosvalue.getDataSize(), result.att.getDataSize());
#pragma omp parallel for #pragma omp parallel for
for (size_t i = 0; i < result.tof.getDataSize(); i++) for (size_t i = 0; i < result.tof.getDataSize(); i++)
{ {
EXPECT_DOUBLE_AE(0,result.tof[i])<<",index:"<<i; EXPECT_DOUBLE_AE(tof[i],result.tof[i])<<",index:"<<i;
EXPECT_DOUBLE_AE(sosvalue[i],result.sosValue[i])<<",index:"<<i; EXPECT_DOUBLE_AE(sosvalue[i],result.sosValue[i])<<",index:"<<i;
EXPECT_TRUE(std::isnan(result.att[i]))<<",index:"<<i; // EXPECT_DOUBLE_AE(result1.att[i],result.att[i])<<",index:"<<i;
// EXPECT_TRUE((std::abs(att[i]-result.att[i])<0.001))<<",att["<<i<<"]"<<att[i]<<", res att["<<i<<"]"<<result.att[1];
} }
} }