diff --git a/src/transmissionReconstruction/detection/detection.cpp b/src/transmissionReconstruction/detection/detection.cpp index 6d60ab6..d8af58f 100644 --- a/src/transmissionReconstruction/detection/detection.cpp +++ b/src/transmissionReconstruction/detection/detection.cpp @@ -84,9 +84,10 @@ namespace Recon { if(gaussWindow) { auto expectedPosWater = (distBlock / expectedSOSWater) * sampleRate + startOffset; + auto windowWidth = calcResult.endSearch-calcResult.startSearch; + #pragma omp parallel for for (size_t i = 0; i < AscanBlock.getDimSize(1); i++) { - auto windowWidth = calcResult.endSearch-calcResult.startSearch; auto t = linspace(-5,5,windowWidth[i]); auto gauss = exp( -0.1 * (transpose(t) ^ 2)); @@ -105,6 +106,7 @@ namespace Recon { } } else{ + #pragma omp parallel for for (size_t i = 0; i < AscanBlock.getDimSize(1); i++) { for (size_t j = calcResult.startSearch[i] - 1; j < calcResult.endSearch[i]; j++) { @@ -118,7 +120,6 @@ namespace Recon { return result; } - //TODO:need test Aurora::Matrix detectAttVectorized(const Aurora::Matrix &Ascan, const Aurora::Matrix &AscanRef, const Aurora::Matrix &distRef, const Aurora::Matrix &sosWaterRef, @@ -170,24 +171,34 @@ namespace Recon { auto sampleRate = aScanReconstructionFrequency; double offsetElectronicSamples = offsetElectronic * sampleRate; Matrix diffStartSearch; + TimeWindowResult timeResult1; + timeResult1.AscanBlockProcessed = AscanBlock; + TimeWindowResult timeResult2; + timeResult2.AscanBlockProcessed = AscanRefBlock; if (useTimeWindowing == 1) { - auto timeResult1 = applyTimeWindowing(AscanBlock, sampleRate, distBlock, sosWaterBlock, expectedSOSWater, offsetElectronicSamples, detectionWindowSOS, minSpeedOfSound, maxSpeedOfSound, gaussWindow); - auto timeResult2= applyTimeWindowing(AscanRefBlock, sampleRate, distBlockRef, sosWaterRefBlock, expectedSOSWater, offsetElectronicSamples, detectionWindowSOS, minSpeedOfSound, maxSpeedOfSound, gaussWindow); + 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; } + 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 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 c1 = ifft(x*conj(y)); - auto c = complex(zeros(mxl+mxl+1,c1.getDimSize(1))); + 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)); + #pragma omp parallel for for (size_t i = 0; i < mxl; i++) { c(i,$) = c1(m2-mxl+i, $); @@ -196,11 +207,12 @@ namespace Recon { 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; + max(c($,i).toMatrix(),All,rowID,colID); + shiftInSamples[i]=-maxlag+colID+rowID; } if (useTimeWindowing) { @@ -208,5 +220,6 @@ namespace Recon { } auto tof = shiftInSamples / sampleRate + distBlock / sosWaterBlock; auto sosValue = distBlock / tof; + return tof; } } diff --git a/test/Detection_Test.cpp b/test/Detection_Test.cpp index 3a362b8..e05b2f7 100644 --- a/test/Detection_Test.cpp +++ b/test/Detection_Test.cpp @@ -73,24 +73,19 @@ TEST_F(Detection_Test, calculateAttenuation) { TEST_F(Detection_Test, applyTimeWindowing) { - MatlabReader m("/home/krad/TestData/timeWindow.mat"); + MatlabReader m("/home/krad/TestData/timeWindow2.mat"); auto AscanBlock = m.read("AscanBlock"); - auto distBlock = m.read("distBlock"); - auto sosBlock = m.read("sosBlock"); - auto AscanBlockProcessed = m.read("AscanBlockProcessed"); + auto distBlock = m.read("dists"); + auto sosBlock = m.read("waterTempBlock"); + auto AscanBlockProcessed = m.read("AscanBlock1"); auto startSearch = m.read("startSearch"); - auto result = Recon::applyTimeWindowing(AscanBlock, 10000000, distBlock, sosBlock, 1482, 0, 0, 1400, 1650, false); + auto result = Recon::applyTimeWindowing(AscanBlock, 10000000, distBlock, sosBlock, 1.511948131508464e+03, 5.2, 1, 1450, 1550, false); + #pragma omp parallel for for (size_t i = 0; i < AscanBlockProcessed.getDataSize(); i++) { EXPECT_DOUBLE_AE(AscanBlockProcessed[i],result.AscanBlockProcessed[i])<<",index:"<