diff --git a/src/main.cxx b/src/main.cxx index 11b47f0..f311e69 100644 --- a/src/main.cxx +++ b/src/main.cxx @@ -9,8 +9,26 @@ #include "common/ceMatchedFilterHandling.h" #include "MatlabReader.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 getLogger(const char* title) { + auto console_sink = std::make_shared(); + console_sink->set_level(spdlog::level::info); + console_sink->set_pattern(fmt::format("[%Y-%m-%d %T .%e][{}] [%^%l%$] %v", title)); + std::shared_ptr logger(new spdlog::logger(title, {console_sink})); + logger->set_level(spdlog::level::info); + logger->flush_on(spdlog::level::info); + return logger; +} + int main() { + auto defaultLogger = getLogger("Main"); + spdlog::set_default_logger(defaultLogger); + SPDLOG_INFO("start"); Recon::startReconstructions(); + SPDLOG_INFO("finish"); return 0; } diff --git a/src/transmissionReconstruction/detection/detection.cpp b/src/transmissionReconstruction/detection/detection.cpp index 329ad0e..ee3e0b8 100644 --- a/src/transmissionReconstruction/detection/detection.cpp +++ b/src/transmissionReconstruction/detection/detection.cpp @@ -151,7 +151,7 @@ namespace Recon { for (size_t i = 0; i < Ascan.getDimSize(1); i++) { 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)); endPosRef[i] = std::ceil(std::min(sizeAscan[1], tof2[i]+detectionWindowATT)); } @@ -194,7 +194,7 @@ namespace Recon { expectedSOSWater, offsetElectronicSamples, detectionWindowSOS, minSpeedOfSound, maxSpeedOfSound, gaussWindow); - diffStartSearch = timeResult1.startSearch - timeResult1.startSearch; + diffStartSearch = timeResult1.startSearch - timeResult2.startSearch; } auto _AscanBlock = timeResult1.AscanBlockProcessed; auto _AscanRefBlock = timeResult2.AscanBlockProcessed; @@ -213,15 +213,17 @@ namespace Recon { auto c1 = real(c_1_2); 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, $); - } - c(mxl + mxl, $) = c1(mxl, $); + 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 + // #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); @@ -248,15 +250,69 @@ namespace Recon { 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; + 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 - 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( - AscanBlock, AscanRefBlock, distRefBlock, sosWaterRefBlock, - result.tof, aScanReconstructionFrequency, offsetElectronic, + _AscanBlock, _AscanRefBlock, distRefBlock, sosWaterRefBlock, + tof, aScanReconstructionFrequency, offsetElectronic, detectionWindowATT); result.tof = tofRel; return result; diff --git a/test/Detection_Test.cpp b/test/Detection_Test.cpp index 2badb7b..b14718f 100644 --- a/test/Detection_Test.cpp +++ b/test/Detection_Test.cpp @@ -63,22 +63,23 @@ TEST_F(Detection_Test, detectTofAndAttMex) { MatlabReader m3("/home/krad/TestData/sosResult.mat"); auto sosvalue = m3.read("sosValue"); + auto tof = m3.read("tofRel"); ASSERT_EQ(sosvalue.getDataSize(), result.tof.getDataSize()); ASSERT_EQ(sosvalue.getDataSize(), result.sosValue.getDataSize()); ASSERT_EQ(sosvalue.getDataSize(), result.att.getDataSize()); #pragma omp parallel for for (size_t i = 0; i < result.tof.getDataSize(); i++) { - EXPECT_DOUBLE_AE(0,result.tof[i])<<",index:"<