diff --git a/src/transmissionReconstruction/dataFilter/dataFilter.cpp b/src/transmissionReconstruction/dataFilter/dataFilter.cpp new file mode 100644 index 0000000..e21880c --- /dev/null +++ b/src/transmissionReconstruction/dataFilter/dataFilter.cpp @@ -0,0 +1,132 @@ +#include "dataFilter.h" + +#include +#include +#include + +#include "Function.h" +#include "Function1D.h" +#include "Function2D.h" +#include "Function3D.h" +#include "Matrix.h" + +using namespace Aurora; +namespace Recon { + Aurora::Matrix calculateSnr(const Aurora::Matrix &aMDataBlock, + double aReferenceNoise) { + auto maxSignal = max(abs(aMDataBlock)); + auto snrBlock = 10 * log(maxSignal / aReferenceNoise, 10); + return snrBlock; + } + + Aurora::Matrix filterTransmissionSensitivityMap( + double sensFilter, const Aurora::Matrix &aVslBlock, + const Aurora::Matrix &aVsnBlock, const Aurora::Matrix &aVrlBlock, + const Aurora::Matrix &aVrnBlock, std::vector &aMSensData) { + if (aMSensData.empty()) { + std::cerr << "aMSensData cant be empty!!!" << std::endl; + return Aurora::Matrix(); + } + auto aMSensData0 = aMSensData[0]; + auto size = zeros(4, 1, 1); + size.getData()[0] = aMSensData0.getDimSize(0); + size.getData()[1] = aMSensData0.getDimSize(1); + size.getData()[2] = aMSensData0.getDimSize(2); + size.getData()[3] = aMSensData.size(); + auto idx = sub2ind(size, {aVrnBlock, aVrlBlock, aVsnBlock, aVslBlock}); + auto transData = zeros(idx.getDataSize(), 1, 1); + for (size_t i = 0; i < idx.getDataSize(); i++) { + auto index = (size_t)(idx.getData()[i] - 1); + auto sliceIndex = index / aMSensData0.getDataSize(); + transData.getData()[i] = + (aMSensData[sliceIndex].getData()[index] >= sensFilter ? 1.0 : 0.0); + } + return transData; + } + + Aurora::Matrix + filterTransmissionAngle(double aAngleLowerLimit, double aAngleUpperLimit, + const Aurora::Matrix &aMSenderNormalBlock, + const Aurora::Matrix &aMReceiverNormalBlock) { + // dot(senderNormalBlock, receiverNormalBlock, + // 1))./(vecnorm(senderNormalBlock, 2, 1) .* vecnorm(receiverNormalBlock, 2, + // 1) + auto transData = ones(1, aMSenderNormalBlock.getDimSize(1)); + auto inbetweenAngle = acosd(dot(aMSenderNormalBlock, aMReceiverNormalBlock) / + (vecnorm(aMSenderNormalBlock, Norm2, 1) * + vecnorm(aMReceiverNormalBlock, Norm2, 1))); + for (size_t i = 0; i < transData.getDataSize(); i++) + { + transData.getData()[i] = (inbetweenAngle.getData()[i] aAngleUpperLimit)?0.0:1.0; + } + return transData; + } + + Aurora::Matrix checkTofDetections(Aurora::Matrix &aVTofValues, const Aurora::Matrix &aVDists, + const Aurora::Matrix &aVSosRef, + double minSpeedOfSound, + double maxSpeedOfSound) + { + auto sosValues = aVDists / (aVTofValues + (aVDists / aVSosRef)); + + auto valid = (sosValues < maxSpeedOfSound) * (sosValues > minSpeedOfSound); + auto minSpeedOfSoundM = minSpeedOfSound+ zeros(sosValues.getDimSize(0),sosValues.getDimSize(1),sosValues.getDimSize(2)); + auto maxSpeedOfSoundM = maxSpeedOfSound+zeros(sosValues.getDimSize(0),sosValues.getDimSize(1),sosValues.getDimSize(2)); + sosValues = max(minSpeedOfSoundM, sosValues); + sosValues = min(maxSpeedOfSoundM, sosValues); + + aVTofValues = (aVDists / sosValues) - (aVDists / aVSosRef); + return sosValues; + } + + Aurora::Matrix filterTransmissionData(int aFilterSensitivity, const Aurora::Matrix &aVslBlock, + const Aurora::Matrix &aVsnBlock, const Aurora::Matrix &aVrlBlock, + const Aurora::Matrix &aVrnBlock, std::vector &aMSensData, + const Aurora::Matrix &aMSenderNormalBlock, + const Aurora::Matrix &aMReceiverNormalBlock, + double* params) + { + switch (aFilterSensitivity) { + case 1: + return filterTransmissionSensitivityMap(params[0], aVslBlock, aVsnBlock, aVrlBlock, aVrnBlock, aMSensData); + case 2: + return filterTransmissionAngle(params[0], params[1], aMSenderNormalBlock, aMReceiverNormalBlock); + } + std::cerr<<"FilterSensitivity value error!"<add){ + valid[i] = 0; + } + } + return valid; + } + +} // namespace Recon + diff --git a/src/transmissionReconstruction/dataFilter/dataFilter.h b/src/transmissionReconstruction/dataFilter/dataFilter.h new file mode 100644 index 0000000..e90ef47 --- /dev/null +++ b/src/transmissionReconstruction/dataFilter/dataFilter.h @@ -0,0 +1,53 @@ +#ifndef __DATAFILTER_H__ +#define __DATAFILTER_H__ +#include "Matrix.h" +#include + +namespace Recon { +Aurora::Matrix calculateSnr(const Aurora::Matrix &aMDataBlock, + double aReferenceNoise); + +Aurora::Matrix filterTransmissionSensitivityMap( + double sensFilter, const Aurora::Matrix &aVslBlock, + const Aurora::Matrix &aVsnBlock, const Aurora::Matrix &aVrlBlock, + const Aurora::Matrix &aVrnBlock, std::vector &aMSensData); + +Aurora::Matrix +filterTransmissionAngle(double aAngleLowerLimit, double aAngleUpperLimit, + const Aurora::Matrix &aMSenderNormalBlock, + const Aurora::Matrix &aMReceiverNormalBlock); + +Aurora::Matrix checkTofDetections(Aurora::Matrix &aVTofValues, + const Aurora::Matrix &aVDists, + const Aurora::Matrix &aVSosRef, + double minSpeedOfSound, + double maxSpeedOfSound); +/** + * filterTransmissionData + * + * @param aFilterSensitivity + * @param aVslBlock + * @param aVsnBlock + * @param aVrlBlock + * @param aVrnBlock + * @param aMSensData + * @param aMSenderNormalBlock + * @param aMReceiverNormalBlock + * @param params + * 如果aFilterSensitivity为1,单个值sensFilter,如果为2,双值angleLowerLimit, + * angleUpperLimit + * @return Aurora::Matrix + */ +Aurora::Matrix filterTransmissionData( + int aFilterSensitivity, const Aurora::Matrix &aVslBlock, + const Aurora::Matrix &aVsnBlock, const Aurora::Matrix &aVrlBlock, + const Aurora::Matrix &aVrnBlock, std::vector &aMSensData, + const Aurora::Matrix &aMSenderNormalBlock, + const Aurora::Matrix &aMReceiverNormalBlock, double *params); + +Aurora::Matrix findDefectTransmissionData(const Aurora::Matrix &aVSNRList,double aSNRDifference); + +//TODO:estimateNoiseValueFromAScans.m + +} // namespace Recon +#endif // __DATAFILTER_H__ \ No newline at end of file diff --git a/test/DataFilter_Test.cpp b/test/DataFilter_Test.cpp new file mode 100644 index 0000000..1ac2b6e --- /dev/null +++ b/test/DataFilter_Test.cpp @@ -0,0 +1,117 @@ +#include +#include + +#include "Function1D.h" +#include "MatlabReader.h" +#include "Matrix.h" +#include "transmissionReconstruction/dataFilter/dataFilter.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 DataFilter_Test : public ::testing::Test { +protected: + static void SetUpDataFilterTester() { + + } + + static void TearDownTestCase() { + } + + void SetUp() { + } + + void TearDown() { + } +}; + +TEST_F(DataFilter_Test, filterTransmissionSensitivityMap) { + double *dataA = new double[4]{1, 1, 1, 1}; + auto slBlock = Aurora::Matrix::fromRawData(dataA, 4, 1); + double *dataB = new double[4]{1, 1, 1, 1}; + auto snBlock = Aurora::Matrix::fromRawData(dataB, 4, 1); + double *data3 = new double[4]{1, 1, 1, 2}; + auto rlBlock = Aurora::Matrix::fromRawData(data3, 4, 1); + double *data4 = new double[4]{1, 2, 3, 1}; + auto rnBlock = Aurora::Matrix::fromRawData(data4, 4, 1); + std::vector a66; + double *data6 = new double[6]{3, 2, 1, 9, 8, 6}; + auto sensData0 = Aurora::Matrix::fromRawData(data6, 3, 2, 1); + a66.push_back(sensData0); + auto result = Recon::filterTransmissionSensitivityMap(0.3, slBlock, snBlock, rlBlock, rnBlock, a66); + EXPECT_EQ(4,result.getDataSize()); + for (size_t i = 0; i < 4; i++) + { + EXPECT_DOUBLE_EQ(1.0,result.getData()[i]); + } + +} + +TEST_F(DataFilter_Test, filterTransmissionAngle) { + double *dataA = new double[12]{0.99,0.99,0.99,0.99,0.10,0.10,0.10,0.10,0,0,0,0}; + auto senderNormalBlock = Aurora::transpose(Aurora::Matrix::fromRawData(dataA, 4, 3)); + double *dataB = new double[12]{0.99,0.99,0.99,0.98,0.10,0.10,0.10,-0.15,0,0,0,0}; + auto receiverNormalBlock = Aurora::transpose(Aurora::Matrix::fromRawData(dataB, 4, 3)); + double angleLowerLimit = 10; + double angleUpperLimit = 180; + + auto result = Recon::filterTransmissionAngle(angleLowerLimit, angleUpperLimit, senderNormalBlock, receiverNormalBlock); + EXPECT_EQ(4,result.getDataSize()); + EXPECT_DOUBLE_EQ(0.0,result.getData()[0]); + EXPECT_DOUBLE_EQ(0.0,result.getData()[1]); + EXPECT_DOUBLE_EQ(0.0,result.getData()[2]); + EXPECT_DOUBLE_EQ(1.0,result.getData()[3]); +} + +TEST_F(DataFilter_Test, checkTofDetections) { + double *dataA = new double[3]{-3.0e-07, -2.6e-06, -2.6e-06}; + auto tofValues = Aurora::Matrix::fromRawData(dataA, 3, 1); + double *dataB = new double[3]{0.238, 0.249, 0.249}; + auto dists = Aurora::Matrix::fromRawData(dataB, 3, 1); + double *data3 = new double[3]{1476.35, 1476.28, 1476.28}; + auto sosRef = Aurora::Matrix::fromRawData(data3, 3, 1); + double minSpeedOfSound = 1400; + double maxSpeedOfSound = 1650; + + auto result = Recon::checkTofDetections(tofValues, dists, sosRef, minSpeedOfSound,maxSpeedOfSound); + EXPECT_EQ(3,result.getDataSize()); + EXPECT_DOUBLE_AE(1479.1025,result.getData()[0]); + EXPECT_DOUBLE_AE(1499.3931,result.getData()[1]); + EXPECT_DOUBLE_AE(1499.3931,result.getData()[2]); +} + +TEST_F(DataFilter_Test, calculateSnr) { + MatlabReader m("/home/krad/TestData/snr.mat"); + + auto snrBlock = m.read("snrBlock"); + auto dataBlock = m.read("dataBlock"); + auto result = Recon::calculateSnr(dataBlock,1.978); + + for (size_t i = 0; i < snrBlock.getDataSize(); i++) + { + EXPECT_DOUBLE_AE(snrBlock.getData()[i],result.getData()[i]); + } +} + +TEST_F(DataFilter_Test, findDefectTransmissionData) { + MatlabReader m("/home/krad/TestData/finddefect.mat"); + double *dataA = new double[3]{1, std::numeric_limits::infinity(), -std::numeric_limits::infinity()}; + auto SNRList = m.read("SNRList"); + auto SNRList2 = Aurora::Matrix::fromRawData(dataA, 3,1,1); + auto result = m.read("valid"); + + auto valid = Recon::findDefectTransmissionData(SNRList2,0.99); + EXPECT_DOUBLE_AE(valid[0],1); + EXPECT_DOUBLE_AE(valid[1],0); + valid = Recon::findDefectTransmissionData(SNRList,0.99); + for (size_t i = 0; i < valid.getDataSize(); i++) + { + EXPECT_DOUBLE_AE(valid[i],0); + } +} + + diff --git a/test/Sensitivity_Test.cpp b/test/Sensitivity_Test.cpp index b98f9af..23f3fc0 100644 --- a/test/Sensitivity_Test.cpp +++ b/test/Sensitivity_Test.cpp @@ -30,12 +30,18 @@ protected: } }; +#include "Function3D.h" + TEST_F(Sensitivity_Test, getSensitivity) { MatlabReader m("/home/krad/TestData/sensitivity.mat"); auto sensmap = m.read("sensmap"); auto senderNormal = m.read("senderNormal"); auto dirToReceiver = m.read("dirVector"); + { auto xyn = Aurora::zeros(1, 3); + xyn.getData()[2] = 1; + xyn = Aurora::repmat(xyn,2304, 1); + } auto output = Recon::getSensitivity(sensmap, senderNormal, dirToReceiver);