diff --git a/src/common/ceMatchedFilterHandling.cpp b/src/common/ceMatchedFilterHandling.cpp new file mode 100644 index 0000000..7645ce2 --- /dev/null +++ b/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; +} + diff --git a/src/common/ceMatchedFilterHandling.h b/src/common/ceMatchedFilterHandling.h new file mode 100644 index 0000000..63d7c6c --- /dev/null +++ b/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