Add ceMatchedFilterHandling in folder common.

src/common/ceMatchedFilterHandling.cpp

@@ -0,0 +1,124 @@
+#include "ceMatchedFilterHandling.h"
+#include "Function.h"
+#include "Function1D.h"
+#include "Function2D.h"
+#include "Function3D.h"
+
+Aurora::Matrix Recon::adaptFrequency(Aurora::Matrix &aCE, double ceSampleFrequency, double requiredFrequency)
+{
+    if (requiredFrequency < ceSampleFrequency){
+        double stride = ceSampleFrequency / (aCE.getDataSize() - 1);
+        double *f2Data = Aurora::malloc(aCE.getDataSize());
+        for (size_t i = 0; i < aCE.getDataSize(); i++)
+        {
+            f2Data[i] = i*stride;
+        }
+        auto f2 = Aurora::Matrix::New(f2Data, aCE.getDataSize(), 1, 1);
+        auto df = fft(aCE);
+        df = df * (f2 < (requiredFrequency / 2));
+        int halfdfLength = (int)std::ceil(((double)aCE.getDataSize()) * 0.5);
+        df.forceReshape(halfdfLength, 1, 1);
+        aCE = Aurora::ifft_symmetric(df, aCE.getDataSize());
+    }
+    double *xData = Aurora::malloc(aCE.getDataSize());
+    for (size_t i = 0; i < aCE.getDataSize(); i++)
+    {
+        xData[i] = (double)(i+1.0);
+    }
+    auto xMatrix = Aurora::Matrix::New(xData,aCE.getDataSize(),1,1); 
+    double stride = ceSampleFrequency /requiredFrequency;
+    double value = 1.0;
+    size_t length = std::ceil(((double)(aCE.getDataSize()-1))/stride);
+    double *x1Data = Aurora::malloc(length);
+    for (size_t i = 0; i < length; i++)
+    {
+        x1Data[i] = value;
+        value+=stride;
+    }
+    auto x1Matrix = Aurora::Matrix::New(x1Data,length,1,1); 
+    
+    aCE = Aurora::interp1(xMatrix,aCE,x1Matrix,Aurora::InterpnMethod::Spline);
+    return aCE;
+}
+
+Aurora::Matrix Recon::preprocessCE(Aurora::Matrix &aCE, double ceSampleFrequency, double requiredFrequency, double expectedLength)
+{
+    adaptFrequency(aCE,ceSampleFrequency,requiredFrequency);
+    aCE = aCE - Aurora::mean(aCE);
+    aCE = aCE / Aurora::max(Aurora::abs(aCE)).getScalar(); 
+    aCE = aCE / Aurora::sum(Aurora::abs(aCE)).getScalar(); 
+
+//     % adapt to AScan length
+// ce(expectedLength) = 0; % padding to size of ascan
+
+// if (length(ce) > expectedLength) % cropping to size of ascan
+//     ce(expectedLength+1:end) = [];
+// end
+
+// ce = ce(:);
+    return aCE;
+}
+
+Aurora::Matrix Recon::createMatchedFilter(const Aurora::Matrix &aCE,
+                                          bool measuredCEused, bool findDefects,
+                                          bool removeOutliersFromCEMeasured, std::string hardwareVersion)
+{
+    Aurora::Matrix mFTime = Aurora::Matrix::copyFromRawData(aCE.getData(), 4000, 2304);
+    if (removeOutliersFromCEMeasured)
+    {
+        auto normSTD = Aurora::std(Aurora::abs(mFTime));
+        Aurora::nantoval(normSTD, 0.0);
+        auto sortSTD = Aurora::sort(normSTD);
+        int t = (int)std::round(0.4 * mFTime.getDimSize(1)) - 1;
+        t = t <= 0 ? 1.0 : t;
+        auto absFTime = abs(mFTime);
+        auto maxAbsFTime = max(absFTime);
+        auto maxFlag = maxAbsFTime < (0.1 * max(maxAbsFTime));
+        auto lessFlag = normSTD < sortSTD(0, t).toMatrix();
+        long maxCol, maxRow;
+        max(normSTD, Aurora::Column, maxRow, maxCol);
+        for (int j = 0; j < mFTime.getDimSize(1); ++j)
+        {
+            if ((bool)(lessFlag.getData()[j]) || (bool)(maxFlag.getData()[j]))
+            {
+                mFTime(Aurora::$, j) = mFTime(Aurora::$, maxCol);
+            }
+        }
+    }
+
+    auto matchedFilter = fft(mFTime);
+    auto sumDiff = Aurora::zeros(1, matchedFilter.getDimSize(1));
+    long minCol = 0, row;
+    {
+        auto absMatchedFilter = Aurora::abs(matchedFilter);
+        auto highNoiseScore = Aurora::mean(absMatchedFilter) * Aurora::std(absMatchedFilter);
+        auto medianNoise = Aurora::median(highNoiseScore);
+        Aurora::min(Aurora::abs(highNoiseScore - Aurora::median(highNoiseScore)), Aurora::Column, row, minCol);
+    }
+    auto maxMatchFilter = matchedFilter(Aurora::$, minCol).toMatrix();
+    for (int k = 0; k < matchedFilter.getDimSize(1); ++k)
+    {
+        sumDiff.getData()[k] = Aurora::sum(Aurora::abs(matchedFilter(Aurora::$, k).toMatrix() - maxMatchFilter)).getData()[0];
+    }
+
+    double meansumDiff = Aurora::mean(sumDiff).getScalar();
+    double stdSumDiff = Aurora::std(sumDiff).getScalar();
+    double sumDiffJ = meansumDiff + 2.596 * stdSumDiff;
+    for (int l = 0; l < sumDiff.getDataSize(); ++l)
+    {
+        if (sumDiff.getData()[l] > sumDiffJ)
+        {
+            matchedFilter(Aurora::$, l) = matchedFilter(Aurora::$, minCol);
+        }
+    }
+
+    if (measuredCEused)
+    {
+        auto mFTime2 = Aurora::real(Aurora::ifft(hardwareVersion == "USCT3Dv3" ? -matchedFilter : matchedFilter));
+        mFTime2 = mFTime2 / Aurora::repmat(Aurora::max(Aurora::abs(mFTime2)), mFTime2.getDimSize(0), 1);
+        mFTime2 = mFTime2 / Aurora::repmat(Aurora::sum(Aurora::abs(mFTime2)), mFTime2.getDimSize(0), 1);
+        matchedFilter = Aurora::fft(mFTime2);
+    }
+    return matchedFilter;
+}
+

src/common/ceMatchedFilterHandling.h

@@ -0,0 +1,11 @@
+#ifndef CE_MATCHED_FILTER_HANDLING_H
+#define CE_MATCHED_FILTER_HANDLING_H
+#include "Matrix.h"
+namespace Recon
+{
+    Aurora::Matrix adaptFrequency(Aurora::Matrix &aCE,double ceSampleFrequency, double requiredFrequency);
+    Aurora::Matrix preprocessCE(Aurora::Matrix &aCE,double ceSampleFrequency, double requiredFrequency, double expectedLength);
+    Aurora::Matrix createMatchedFilter(const Aurora::Matrix &aCE,
+        bool measuredCEused, bool findDefects, bool removeOutliersFromCEMeasured, std::string hardwareVersion);
+}
+#endif // !CE_MATCHED_FILTER_HANDLING_H