From 9037838b5e35f5f604215e6d95216ac0038508b8 Mon Sep 17 00:00:00 2001
From: kradchen <chenjuijun@realcan.cn>
Date: Fri, 28 Apr 2023 13:22:53 +0800
Subject: [PATCH] CreateMatchFilter in main

---
 src/main.cxx | 96 +++++++++++++++++++++++++++++++++-------------------
 1 file changed, 62 insertions(+), 34 deletions(-)

diff --git a/src/main.cxx b/src/main.cxx
index d240f0e..93f2a88 100644
--- a/src/main.cxx
+++ b/src/main.cxx
@@ -6,7 +6,6 @@
 
 #include <complex>
 
-
 #include "Matrix.h"
 #include "Function.h"
 #include "Function1D.h"
@@ -14,63 +13,92 @@
 #include "Function3D.h"
 #include "MatlabReader.h"
 
-int main() {
+int main()
+{
     MatlabReader m;
     Input i;
     MatchedFilter o;
 
-    m.read(&i,&o,"/home/krad/TestData/testCreateMatchedFilter.mat");
-    bool measuredCEused = true,findDefects = i.mParams->mFindDefects,removeOutliers = i.mRemoveOutliersFromCEMeasured;
+    m.read(&i, &o, "/home/krad/TestData/testCreateMatchedFilter.mat");
+    bool measuredCEused = true, findDefects = i.mParams->mFindDefects, removeOutliers = i.mRemoveOutliersFromCEMeasured;
 
-    Aurora::Matrix mFTime = Aurora::Matrix::fromRawData(i.mCe,4000,2304);
-    if (removeOutliers){
+    Aurora::Matrix mFTime = Aurora::Matrix::fromRawData(i.mCe, 4000, 2304);
+    if (removeOutliers)
+    {
         auto normSTD = Aurora::std(Aurora::abs(mFTime));
-        Aurora::nantoval(normSTD,0.0);
+        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;
+        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 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 = 998, row;
+
     auto absMatchedFilter = abs(matchedFilter);
-//    auto highNoiseScore = mean(absMatchedFilter) * Aurora::std(absMatchedFilter);
-    auto highNoiseScore = mean( abs(matchedFilter)) * Aurora::std(absMatchedFilter);
-    printf("highNoiseScore[998]: %f4, highNoiseScore[2053]:%f4\r\n",highNoiseScore.getData()[998],highNoiseScore.getData()[2053]);
+    auto highNoiseScore = mean(absMatchedFilter) * Aurora::std(absMatchedFilter);
+    // auto highNoiseScore = mean( abs(matchedFilter)) * Aurora::std(absMatchedFilter);
+    printf("highNoiseScore[998]: %f4, highNoiseScore[2053]:%f4\r\n", highNoiseScore.getData()[998], highNoiseScore.getData()[2053]);
     auto medianNoise = median(highNoiseScore);
-    printf("median: %f4\r\n",medianNoise.getScalar());
-    long minCol,row;
+    printf("median: %f , should be: 1724817516.8468074\r\n", medianNoise.getScalar());
+
     min(abs(highNoiseScore - median(highNoiseScore)), Aurora::Column, row, minCol);
     //
-    minCol = 998;
-    auto sumDiff = Aurora::zeros(1, matchedFilter.getDimSize(1));
+    // minCol = 998;
+
     auto maxMatchFilter = matchedFilter(Aurora::$, minCol).toMatrix();
-    for (int k = 0; k < matchedFilter.getDimSize(1); ++k) {
-        sumDiff.getData()[k] = sum(abs(matchedFilter(Aurora::$, k).toMatrix()- maxMatchFilter)).getData()[0];
+    for (int k = 0; k < matchedFilter.getDimSize(1); ++k)
+    {
+        sumDiff.getData()[k] = sum(abs(matchedFilter(Aurora::$, k).toMatrix() - maxMatchFilter)).getData()[0];
     }
 
-    auto indexe = sumDiff > (mean(sumDiff)+2.596 * Aurora::std(sumDiff).getScalar());
-    for (int l = 0; l < indexe.getDataSize(); ++l) {
-        if ((bool)indexe.getData()[l]){
-            matchedFilter(Aurora::$,l) = matchedFilter(Aurora::$,minCol);
+    // printf("meanSumDiff\r\n");
+    // auto meanSumDiff = mean(sumDiff);
+    // printf("meanSumDiff finish\r\n");
+    {
+        double meansumDiff = mean(sumDiff).getScalar() ;
+        double stdSumDiff =  Aurora::std(sumDiff).getScalar();
+        double sumDiffJ =meansumDiff + 2.596 * stdSumDiff;
+        // auto indexe = sumDiff > sumDiffJ;
+        printf("indexe finish\r\n");
+
+        for (int l = 0; l < sumDiff.getDataSize(); ++l)
+        {
+            
+            if (sumDiff.getData()[l]> sumDiffJ)
+            {
+                matchedFilter(Aurora::$, l) = matchedFilter(Aurora::$, minCol);
+            }
         }
+        printf("\r\n");
     }
-    if (measuredCEused){
-        auto mFTime2 = ifft(-matchedFilter);
-        max(Aurora::abs(mFTime2)).printfShape();
-        mFTime2 = mFTime2/repmat(max(Aurora::abs(mFTime2)),mFTime2.getDimSize(0),1);
-        mFTime2 = mFTime2/repmat(sum(Aurora::abs(mFTime2)),mFTime2.getDimSize(0),1);
+    
+
+    if (measuredCEused)
+    {
+        auto mFTime2 = real(ifft(-matchedFilter));
+        mFTime2 = mFTime2 / repmat(max(Aurora::abs(mFTime2)), mFTime2.getDimSize(0), 1);
+        mFTime2 = mFTime2 / repmat(sum(Aurora::abs(mFTime2)), mFTime2.getDimSize(0), 1);
         matchedFilter = fft(mFTime2);
     }
-    printf("run end\r\n");
+    for (size_t i = 0; i < 10000; i += 500)
+    {
+        printf("index :%d,origin output:%f4,%f4, output:%f4,%f4\r\n",
+               i, o.mReal[i], o.mImag[i], matchedFilter.getData()[i * 2], matchedFilter.getData()[i * 2 + 1]);
+    }
+
     return 0;
 }