diff --git a/src/transmissionReconstruction/detection/detection.cpp b/src/transmissionReconstruction/detection/detection.cpp new file mode 100644 index 0000000..1d850bd --- /dev/null +++ b/src/transmissionReconstruction/detection/detection.cpp @@ -0,0 +1,210 @@ +#include "detection.h" + +#include "Function1D.h" +#include "Function2D.h" +#include "Function3D.h" +#include "Matrix.h" + +using namespace Aurora; +namespace Recon { + Matrix analyzeDetections(const Matrix &slBlockTotal, const Matrix &snBlockTotal, const Matrix &rlBlockTotal, const Matrix &rnBlockTotal, const Matrix &tofDataTotal){ + //TODO: 待完成,暂时用不到 + return Matrix(); + } + + + Matrix calculateSOSOffset(const Matrix &sosBlock, double referenceSOS, const Matrix &distBlock, double sampleRate){ + auto offset = (distBlock / sosBlock - distBlock / referenceSOS) * sampleRate; + return offset; + } + + + Matrix calculateAttenuation(const Matrix &ascans, const Matrix &startPos, const Matrix &endPos, const Matrix &ascansRef, const Matrix &startPosRef, const Matrix &endPosRef) + { + auto ascans2 = zeros(ascans.getDimSize(0), ascans.getDimSize(1)); + auto ascansRef2 = zeros(ascansRef.getDimSize(0), ascansRef.getDimSize(1)); + #pragma omp parallel for + for (size_t i = 0; i < ascans.getDimSize(1); i++) + { + for (size_t j = startPos[i]-1; j < endPos[i]; j++) + { + ascans2(j, i) = ascans(j, i); + } + for (size_t j = startPosRef[i]-1; j < endPosRef[i]; j++) + { + ascansRef2(j, i) = ascansRef(j, i); + } + } + + auto pulseEnergy = sum(ascans2^2); + auto pulseEnergyEmpty = sum(ascansRef2^2); + + return log(pulseEnergyEmpty/pulseEnergy); + } + + SearchPosition calculateStarEndSearchPosition(const Matrix &aVDistBlock, + double minSpeedOfSound, double maxSpeedOfSound, + double sampleRate, double maxSample, + const Matrix &aVSosOffsetBlock, + double startOffset, double segmentLenOffset) + { + auto startSearch = floor((aVDistBlock / maxSpeedOfSound)*sampleRate+aVSosOffsetBlock+startOffset); + for (size_t i = 0; i < startSearch.getDataSize(); i++) + { + startSearch[i] = (uint)(startSearch[i]); + if(startSearch[i]<1) startSearch[i]=1; + else if(startSearch[i]>maxSample) startSearch[i] = maxSample; + } + + auto endSearch = ceil(aVDistBlock/minSpeedOfSound*sampleRate+aVSosOffsetBlock+startOffset+segmentLenOffset); + for (size_t i = 0; i < endSearch.getDataSize(); i++) + { + endSearch[i] = (uint)(endSearch[i]); + if(endSearch[i]<1) endSearch[i]=1; + else if(endSearch[i]>maxSample) endSearch[i] = maxSample; + } + SearchPosition result; + result.startSearch = startSearch; + result.endSearch = endSearch; + return result; + } + + TimeWindowResult applyTimeWindowing(const Aurora::Matrix &AscanBlock, double sampleRate, + const Aurora::Matrix &distBlock, const Aurora::Matrix &sosBlock, + double expectedSOSWater, double startOffset, double segmentLenOffset, + double minSpeedOfSound, double maxSpeedOfSound, bool gaussWindow) + { + auto sosOffset = calculateSOSOffset(sosBlock, expectedSOSWater, distBlock, sampleRate); + + auto calcResult = calculateStarEndSearchPosition(distBlock, minSpeedOfSound, maxSpeedOfSound, sampleRate,AscanBlock.getDimSize(0), sosOffset, startOffset, segmentLenOffset); + + auto AscanBlockProcessed = zeros(AscanBlock.getDimSize(0),AscanBlock.getDimSize(1)); + + if(gaussWindow) + { + auto expectedPosWater = (distBlock / expectedSOSWater) * sampleRate + startOffset; + 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)); + + for (size_t j = calcResult.startSearch[i]-1; j < calcResult.endSearch[i]; j++) + { + AscanBlockProcessed(j, i) = AscanBlock(j, i); + } + auto temp = zeros(AscanBlockProcessed.getDimSize(0), 1); + size_t end = round(expectedPosWater[i])-floor(windowWidth[i]/2)+gauss.getDataSize()-1; + for (size_t k = floor(windowWidth[i]/2); k < end; k++) + { + temp(round(expectedPosWater[i])-k) = gauss; + } + AscanBlockProcessed($,i) = AscanBlockProcessed($,i).toMatrix() * temp; + } + } + else{ + for (size_t i = 0; i < AscanBlock.getDimSize(1); i++) { + for (size_t j = calcResult.startSearch[i] - 1; + j < calcResult.endSearch[i]; j++) { + AscanBlockProcessed(j, i) = AscanBlock(j, i); + } + } + } + TimeWindowResult result; + result.startSearch = calcResult.startSearch; + result.AscanBlockProcessed = AscanBlockProcessed; + return result; + } + + //TODO:need test + 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) + { + auto sizeAscan = size(Ascan); + + auto sampleRate = aScanReconstructionFrequency; + double offsetElectronicSamples = offsetElectronic * sampleRate; + + auto envelopeOfAScan = abs(hilbert(Ascan)); + auto envelopeOfReferenceAScan = abs(hilbert(AscanRef)); + auto tof2 = (distRef / sosWaterRef) * sampleRate + offsetElectronicSamples; + + auto startPos = zeros(Ascan.getDimSize(1), 1); + auto endPos = zeros(Ascan.getDimSize(1), 1); + auto startPosRef = zeros(Ascan.getDimSize(1), 1); + 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)); + } + return calculateAttenuation(envelopeOfAScan,startPos,endPos,envelopeOfReferenceAScan,startPosRef,endPosRef); + } + + int nextpow2(unsigned int value){ + unsigned int compareValue = 2; + int result = 1; + while(compareValue<=value){ + compareValue=compareValue<<1; + result++; + } + return result; + } + //TODO:need test + Aurora::Matrix 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) + { + auto sampleRate = aScanReconstructionFrequency; + double offsetElectronicSamples = offsetElectronic * sampleRate; + Matrix diffStartSearch; + 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); + + diffStartSearch = timeResult1.startSearch - timeResult1.startSearch; + } + + 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 = zeros(mxl+mxl+1,c1.getDimSize(1)); + for (size_t i = 0; i < mxl; i++) + { + c(i,$) = c1(m2-mxl+(1:mxl), $); + c(i+mxl,$) = c1(i, $); + } + c(mxl+mxl,$) = c1(mxl, $); + + auto shiftInSamples =zeros(1,c1.getDimSize(1)); + 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; + } + if (useTimeWindowing) + { + shiftInSamples = shiftInSamples - diffStartSearch; + } + auto tof = shiftInSamples / sampleRate + distBlock / sosWaterBlock; + auto sosValue = distBlock / tof; + } +} diff --git a/src/transmissionReconstruction/detection/detection.h b/src/transmissionReconstruction/detection/detection.h new file mode 100644 index 0000000..158ef58 --- /dev/null +++ b/src/transmissionReconstruction/detection/detection.h @@ -0,0 +1,47 @@ +#ifndef __TRANS_DETECTION_H__ +#define __TRANS_DETECTION_H__ +#include "Matrix.h" +namespace Recon { +struct SearchPosition { + Aurora::Matrix startSearch; + Aurora::Matrix endSearch; +}; +struct TimeWindowResult { + Aurora::Matrix startSearch; + Aurora::Matrix AscanBlockProcessed; +}; + +struct DetectTofResult { + Aurora::Matrix tof; + Aurora::Matrix sosValue; +}; +Aurora::Matrix calculateAttenuation(const Aurora::Matrix &ascans, + const Aurora::Matrix &startPos, + const Aurora::Matrix &endPos, + const Aurora::Matrix &ascansRef, + const Aurora::Matrix &startPosRef, + const Aurora::Matrix &endPosRef); + +SearchPosition +calculateStarEndSearchPosition(const Aurora::Matrix &aVDistBlock, + double minSpeedOfSound, double maxSpeedOfSound, + double sampleRate, double maxSample, + const Aurora::Matrix &aVSosOffsetBlock, + double startOffset, double segmentLenOffset); + +TimeWindowResult applyTimeWindowing( + const Aurora::Matrix &AscanBlock, double sampleRate, + const Aurora::Matrix &distBlock, const Aurora::Matrix &sosBlock, + double expectedSOSWater, double startOffset, double segmentLenOffset, + double minSpeedOfSound, double maxSpeedOfSound, bool gaussWindow); + +// 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); + +} // namespace Recon + +#endif // __DETECTION_H__ \ No newline at end of file diff --git a/test/Detection_Test.cpp b/test/Detection_Test.cpp new file mode 100644 index 0000000..75722dc --- /dev/null +++ b/test/Detection_Test.cpp @@ -0,0 +1,85 @@ +#include +#include + +#include "Function1D.h" +#include "MatlabReader.h" +#include "Matrix.h" +#include "transmissionReconstruction/detection/detection.h" + + + + +inline double fourDecimalRound(double src){ + return round(src*10000.0)/10000.0; +} + +#define EXPECT_DOUBLE_AE(valueA,valueB)\ + EXPECT_DOUBLE_EQ(fourDecimalRound(valueA),fourDecimalRound(valueB)) + +class Detection_Test : public ::testing::Test { +protected: + static void SetUpDetectionTester() { + + } + + static void TearDownTestCase() { + } + + void SetUp() { + } + + void TearDown() { + } +}; + +TEST_F(Detection_Test, calculateStarEndSearchPosition) { + + auto distBlock = Aurora::Matrix::fromRawData(new double[3]{0.22, 0.21, 0.11}, 3, 1); + auto sosOffsetBlock = Aurora::Matrix::fromRawData(new double[3]{-0.8, 0, 0.9}, 3, 1); + + auto result = Recon::calculateStarEndSearchPosition(distBlock, 1400.0, 1650.0, 10000000, 9999, sosOffsetBlock,97.3,250); + EXPECT_EQ(3,result.endSearch.getDataSize()); + EXPECT_EQ(3,result.startSearch.getDataSize()); + EXPECT_DOUBLE_AE(1429,result.startSearch[0]); + EXPECT_DOUBLE_AE(1370,result.startSearch[1]); + EXPECT_DOUBLE_AE(764,result.startSearch[2]); + EXPECT_DOUBLE_AE(1918,result.endSearch[0]); + EXPECT_DOUBLE_AE(1848,result.endSearch[1]); + EXPECT_DOUBLE_AE(1134,result.endSearch[2]); + +} + +TEST_F(Detection_Test, calculateAttenuation) { + + MatlabReader m("/home/krad/TestData/calcAtt.mat"); + + auto ascans = m.read("ascans"); + auto ascansRef = m.read("ascansRef"); + auto endPos = m.read("endPos"); + auto endPosRef = m.read("endPosRef"); + auto startPos = m.read("startPos"); + auto startPosRef = m.read("startPosRef"); + auto att = m.read("att"); + auto result = Recon::calculateAttenuation(ascans, startPos, endPos, ascansRef, startPosRef, endPosRef); + for (size_t i = 0; i < att.getDataSize(); i++) + { + EXPECT_DOUBLE_AE(att[i],result[i]); + } +} + +TEST_F(Detection_Test, applyTimeWindowing) { + + MatlabReader m("/home/krad/TestData/timeWindow.mat"); + + auto AscanBlock = m.read("AscanBlock"); + auto distBlock = m.read("distBlock"); + auto sosBlock = m.read("sosBlock"); + auto AscanBlockProcessed = m.read("AscanBlockProcessed"); + auto startSearch = m.read("startSearch"); + auto result = Recon::applyTimeWindowing(AscanBlock, 10000000, distBlock, sosBlock, 1482, 0, 0, 1400, 1650, false); + for (size_t i = 0; i < AscanBlockProcessed.getDataSize(); i++) + { + EXPECT_DOUBLE_AE(AscanBlockProcessed[i],result.AscanBlockProcessed[i])<<",index:"<