Add TransmissionReconstruction Datafilter.

2023-05-16 13:18:59 +08:00
parent 21750be9be
commit a0033d9f8a
4 changed files with 308 additions and 0 deletions
--- a/src/transmissionReconstruction/dataFilter/dataFilter.cpp
+++ b/src/transmissionReconstruction/dataFilter/dataFilter.cpp
@@ -0,0 +1,132 @@
+#include "dataFilter.h"
+
+#include <cmath>
+#include <cstddef>
+#include <iostream>
+
+#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<Aurora::Matrix> &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]<aAngleLowerLimit || 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<Aurora::Matrix> &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!"<<std::endl;
+        return Aurora::Matrix();
+    }
+
+    Aurora::Matrix findDefectTransmissionData(const Aurora::Matrix &aVSNRList,double aSNRDifference)
+    {
+        auto valid = auroraNot(Aurora::isnan(aVSNRList)); 
+        double meanSNR = 0.0;
+        
+        auto finite_SNRList = Aurora::isfinite(aVSNRList);
+        int count = (int)sum(finite_SNRList).getScalar();
+        double* std_SNRListData = Aurora::malloc(count);
+        int j = 0;
+        for (size_t i = 0; i < aVSNRList.getDataSize(); i++)
+        {
+            if (finite_SNRList.getData()[i] == 1.0){
+                meanSNR+=aVSNRList[i];
+                std_SNRListData[j++] = aVSNRList[i];
+            }
+        }
+        Aurora::Matrix std_SNRList = Aurora::Matrix::New(std_SNRListData,count,1,1);
+        std_SNRList = Aurora::std(std_SNRList);
+        double localSNRDifference  = 2 * std_SNRList.getScalar();
+        aSNRDifference = localSNRDifference < aSNRDifference?localSNRDifference:aSNRDifference;
+        double sub = meanSNR - aSNRDifference;
+        double add = meanSNR + aSNRDifference;
+        for (size_t i = 0; i < aVSNRList.getDataSize(); i++)
+        {
+            double value = aVSNRList[i];
+            if (value<sub || value>add){
+                valid[i] = 0;
+            }
+        }
+        return valid;
+    }
+
+} // namespace Recon
+
--- a/src/transmissionReconstruction/dataFilter/dataFilter.h
+++ b/src/transmissionReconstruction/dataFilter/dataFilter.h
@@ -0,0 +1,53 @@
+#ifndef __DATAFILTER_H__
+#define __DATAFILTER_H__
+#include "Matrix.h"
+#include <vector>
+
+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<Aurora::Matrix> &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<Aurora::Matrix> &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__
--- a/test/DataFilter_Test.cpp
+++ b/test/DataFilter_Test.cpp
@@ -0,0 +1,117 @@
+#include <gtest/gtest.h>
+#include <limits>
+
+#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<Aurora::Matrix> 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<double>::infinity(), -std::numeric_limits<double>::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);
+    }
+}
+
+
--- 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);