From 48f3d4d7f83c15a6d4d8f22b1eecec2a1978077b Mon Sep 17 00:00:00 2001
From: kradchen <chenjuijun@realcan.cn>
Date: Tue, 13 Jun 2023 14:24:05 +0800
Subject: [PATCH] preprocessAScanBlockForReflection Test passed (left with
 Accuracy problem)

---
 .../preprocessAScanBlockForReflection.cpp     | 439 ++++++++++++++++++
 .../preprocessAScanBlockForReflection.h       |  26 ++
 test/Reconstruction_Test.cpp                  |  49 +-
 3 files changed, 506 insertions(+), 8 deletions(-)
 create mode 100644 src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp
 create mode 100644 src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.h

diff --git a/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp b/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp
new file mode 100644
index 0000000..f7dfc7d
--- /dev/null
+++ b/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp
@@ -0,0 +1,439 @@
+#include "preprocessAScanBlockForReflection.h"
+
+#include "Function.h"
+#include "Function1D.h"
+#include "Function2D.h"
+#include "Function3D.h"
+#include "Matrix.h"
+#include "common/dataBlockCreation/removeDataFromArrays.h"
+#include "config/config.h"
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <cstring>
+#include <iostream>
+#include <mkl_cblas.h>
+#include <mkl_vml_functions.h>
+#include <vector>
+namespace Recon {
+
+    Aurora::Matrix _createDiffVector(const Aurora::Matrix & aMatrix, double beginValue,double endValue){
+        auto r = Aurora::zeros(1,aMatrix.getDataSize()+1);
+        r[0] = beginValue;
+        for (size_t i = 1; i < aMatrix.getDataSize(); i++)
+        {   
+            r[i] = aMatrix[i] - aMatrix[i-1];
+        }
+        r[aMatrix.getDataSize()] = endValue;
+        return r;
+    }
+
+    /**
+     * @brief 
+     * 
+     * @param aMatrix 
+     * @param op add 0, sub 1
+     * @return Aurora::Matrix 
+     */
+    Aurora::Matrix _blockOp(const Aurora::Matrix & aMatrix,int op){
+        if (aMatrix.isComplex()){
+            std::cerr<<"MatrixColDiff not support complex!"<<std::endl;
+            return Aurora::Matrix();
+        }
+        else{
+            if (aMatrix.isVector())
+            {
+                Aurora::Matrix result = Aurora::zeros(aMatrix.getDimSize(0)>1?aMatrix.getDataSize()-1:1,
+                            aMatrix.getDimSize(0)>1?1:aMatrix.getDimSize(1)-1);
+                if (op == 1)vdSubI(result.getDataSize(),aMatrix.getData() , 1, aMatrix.getData()+1,1,result.getData(),1);
+                if (op == 0)vdAddI(result.getDataSize(),aMatrix.getData()+1 , 1, aMatrix.getData(),1,result.getData(),1);
+                return result;
+            }
+            else{
+                std::cerr<<"un expect operation!"<<std::endl;
+                return Aurora::Matrix();
+                // Aurora::Matrix result = Aurora::zeros(aMatrix.getDimSize(0)-1,aMatrix.getDimSize(1));
+                // #pragma omp parallel for
+                // for (size_t i = 0; i < aMatrix.getDimSize(1); i++)
+                // {
+                //     if (op == 1){
+                //         vdSubI(result.getDimSize(0),
+                //             aMatrix.getData()+i*aMatrix.getDimSize(0)+1, 1,
+                //             aMatrix.getData()+i*aMatrix.getDimSize(0)+1, 1,
+                //             result.getData()+i*aMatrix.getDimSize(0),1);
+                //     }
+                //     if (op == 0){
+                //         vdAddI(result.getDimSize(0),
+                //             aMatrix.getData()+i*aMatrix.getDimSize(0)+1, 1,
+                //             aMatrix.getData()+i*aMatrix.getDimSize(0)+1, 1,
+                //             result.getData()+i*aMatrix.getDimSize(0),1);
+                //     }
+
+            }
+        }
+
+    }
+
+    Aurora::Matrix generateGaussWindow(size_t winLength)
+    {
+        auto n = Aurora::linspace(0,5,winLength);
+        return exp( -0.2 * (transpose(n)^ 2));
+    }
+
+
+
+    //TODO:untested
+    Aurora::Matrix offsetSignalFourier(const PreComputes &preComputes, int nfft)
+    {
+         
+        double dt = -(preComputes.offset/preComputes.timeInterval);
+        double binStart = floor((double)nfft/2);
+        auto temp = Aurora::zeros(1,nfft);
+        for (size_t i = 0; i < temp.getDataSize(); i++)
+        {
+            temp[i]=i;
+        }
+        temp= temp-binStart;
+        temp = Aurora::transpose(temp);
+        Aurora::ifftshift(temp);
+        auto fftBin = -1*dt*2*M_PI*temp/nfft;
+        auto negCfftbin = Aurora::complex(Aurora::zeros(fftBin.getDimSize(0),fftBin.getDimSize(1),fftBin.getDimSize(2)));
+        cblas_dcopy(fftBin.getDataSize(), fftBin.getData(), 1, negCfftbin.getData()+1, 2);
+        auto exponent = exp(negCfftbin);
+        if (!exponent.isComplex())exponent = complex(exponent);
+        return exponent;
+    }
+
+
+    std::vector<Aurora::Matrix> performSignalProcessing(
+    const Aurora::Matrix &blockedAScans, 
+    const Aurora::Matrix &blockedSL, 
+    const Aurora::Matrix &blockedRL, 
+    const Aurora::Matrix &blockedSenderPosition,
+    const Aurora::Matrix &blockedReceiverPosition,
+    const Aurora::Matrix &blockedGain, const Aurora::Matrix &blockedChannels,
+    ExpInfo &info, const PreComputes &preComputes) {
+        std::vector<Aurora::Matrix> result;
+        auto t = size(blockedAScans);
+        size_t nSamples = (int)t[0];
+        size_t nAScans = (int)t[1];
+        auto valid = Aurora::ones(blockedAScans.getDimSize(1),1);
+        {
+            auto blockedGain_t = blockedGain.deepCopy();
+            Aurora::nantoval(blockedGain_t, 0);
+            valid = valid*blockedGain_t;
+        }
+        
+        auto _blockedAScans = blockedAScans/blockedGain;
+        if (reflectParams::removeDCOffset == 1)
+	    {
+            _blockedAScans = _blockedAScans-Aurora::mean(_blockedAScans);
+        }
+        int winLength = 100; 
+        if (reflectParams::removeDCOffset == 1)
+        {
+            auto meanAS = mean(_blockedAScans.block(0,winLength-1,_blockedAScans.getDimSize(0) -1 -winLength));
+            _blockedAScans = _blockedAScans-meanAS.getScalar();
+        }
+        _blockedAScans.setBlockValue(0,0,winLength-1,0);
+        _blockedAScans.setBlockValue(0,_blockedAScans.getDimSize(0)-winLength-1,_blockedAScans.getDimSize(0)-1,0);
+        auto maxVal = max(abs(_blockedAScans));
+        auto stdVal = Aurora::std(_blockedAScans);
+        auto meanVal = Aurora::mean(abs(_blockedAScans));
+        auto ascanMapValue = meanVal*stdVal;
+        auto temp = (stdVal == 0) * (maxVal != 0);
+        Aurora::compareSet(valid,temp,1,0,Aurora::EQ);
+        if (reflectParams::findDefects==1)
+        {
+            valid = valid*(maxVal/meanVal >= reflectParams::epsilon);
+        } 
+        // debugOutput.maxVal = maxVal;
+        // debugOutput.stdVal = stdVal;
+        // debugOutput.meanVal = meanVal;
+        // debugOutput.snr = maxVal ./ meanVal;
+        // debugOutput.snrPass =  (maxVal./meanVal) < flags.dataSelection.epsilon;
+        if (reflectParams::suppressSameHead == 1)
+        {
+            for (size_t i = 0; i < blockedSL.getDataSize(); i++)
+            {
+                if (blockedSL[i] == blockedRL[i])
+                {
+                    auto begin_ptr = _blockedAScans.getData()+i*_blockedAScans.getDimSize(0);
+                    int count = reflectParams::suppressSameHeadLength* sizeof(double);
+                    std::memset(begin_ptr,0,count);
+                }
+            }
+
+        }
+        auto fx = fft(_blockedAScans);
+        Aurora::Matrix fh;
+        if (reflectParams::useCorrelation && reflectParams::matchedFilterCeAScan)
+        {
+            if (!preComputes.measuredCEUsed)
+            {
+                fh = preComputes.matchedFilter.block(1, 0, nAScans-1);
+            }
+            else{
+                fh = preComputes.matchedFilter(Aurora::$,blockedChannels[0]-1).toMatrix();
+            }
+            double* value1 = Aurora::malloc(fx.getDataSize());
+            vdMulI(fx.getDataSize(), fx.getData(), 2, fh.getData(), 2, value1, 1);
+            double* value2 = Aurora::malloc(fx.getDataSize());
+            vdMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData() + 1, 2, value2, 1);
+            double* realData = Aurora::malloc(fx.getDataSize());
+            vdAdd(fx.getDataSize(), value1, value2, realData);
+            Aurora::Matrix real = Aurora::Matrix::New(realData, fx.getDimSize(0), fx.getDimSize(1));
+
+            vdMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData(), 2, value1, 1);
+            vdMulI(fx.getDataSize(), fx.getData(), 2, fh.getData() + 1, 2, value2, 1);
+            double* imagData = Aurora::malloc(fx.getDataSize());
+            vdSub(fx.getDataSize(), value1, value2, imagData);
+            Aurora::Matrix image = Aurora::Matrix::New(imagData, fx.getDimSize(0), fx.getDimSize(1));
+
+            double* complexData = Aurora::malloc(real.getDataSize(), true);
+            cblas_dcopy(real.getDataSize(), real.getData(), 1 , complexData ,2);
+            cblas_dcopy(image.getDataSize(), image.getData(), 1 , complexData + 1 ,2);
+            Aurora::Matrix complex = Aurora::Matrix::New(complexData, real.getDimSize(0), real.getDimSize(1), 1, Aurora::Complex);
+            _blockedAScans = Aurora::real(ifft(complex));
+
+            _blockedAScans.setBlockValue(0, 0, winLength-1, 0);
+            _blockedAScans.setBlockValue(0, _blockedAScans.getDimSize(0)-winLength, _blockedAScans.getDimSize(0)-1, 0);
+            _blockedAScans = fft(_blockedAScans);
+        }
+        else{
+            _blockedAScans = fx;
+        }
+        auto exponent = offsetSignalFourier(preComputes, nSamples);
+        // exponent.forceReshape(1, exponent.getDataSize(), 1);
+        _blockedAScans =_blockedAScans*exponent;
+        _blockedAScans = real(ifft(_blockedAScans));
+        if (reflectParams::removeTransmissionSignal == 1) {
+            auto distanceEPosRPos =
+                sqrt(
+                    ((blockedSenderPosition(0, Aurora::$).toMatrix() - blockedReceiverPosition(0, Aurora::$).toMatrix())^2) + 
+                      ((blockedSenderPosition(1, Aurora::$).toMatrix() - blockedReceiverPosition(1, Aurora::$).toMatrix())^2) +
+                       ((blockedSenderPosition(2, Aurora::$).toMatrix() -blockedReceiverPosition(2, Aurora::$).toMatrix())^2));
+            auto windowLength = reflectParams::windowLength;
+            auto start = round((distanceEPosRPos/reflectParams::expectedUSSpeedRange[1])*reflectParams::aScanReconstructionFrequency);
+            auto ende  = round((distanceEPosRPos/reflectParams::expectedUSSpeedRange[0])*reflectParams::aScanReconstructionFrequency);
+            auto area    = ceil(reflectParams::expectedPulseLength*info.Wavelength+1);
+            start   = start- round(area*0.1);
+            ende    = ende + round(area*0.9);
+            Aurora::compareSet(start, 1, 2, Aurora::NG);
+            Aurora::compareSet(start,nSamples,nSamples,Aurora::GT);
+            Aurora::compareSet(ende,nSamples,nSamples,Aurora::GT);
+            
+            for (size_t idx = 0; idx < nAScans; idx++)
+            {
+                auto partData = _blockedAScans(Aurora::$,idx).toMatrix().block(0, start[idx]-1, ende[idx]-1);
+        
+                auto l_partData = partData.getDataSize();
+                auto energyMove = _createDiffVector(partData,partData[0],0);
+                energyMove     = _blockOp(energyMove,0);
+                energyMove     = (energyMove^2)*0.5;
+                auto energyPot = partData^2;
+        
+                auto energySum = energyMove + energyPot;
+        
+                auto mean_energySum = Aurora::zeros(l_partData + windowLength,1);
+                for (size_t i = 0; i < windowLength; i++)
+                {
+                    double* begin = mean_energySum.getData() + i;
+                    int length = l_partData;
+                    vdAddI(length, energySum.getData(), 1, begin, 1, begin, 1);
+                }
+                mean_energySum = mean_energySum.block(0,windowLength/2-1,mean_energySum.getDataSize()-1-round((windowLength-0.001)/2));
+                
+                auto g = _createDiffVector(mean_energySum,0,0);
+                g = _blockOp(g,0);
+                auto f = _createDiffVector(g,g[0],0);
+                f = _blockOp(f,0);
+                Aurora::compareSet(f, 0, 0, Aurora::LT);
+                Aurora::compareSet(f,g,0,0,Aurora::NG);
+                Aurora::compareSet(f,0,0,Aurora::NG);
+                auto mean_energySum_f = f;
+                long scheitelRow=0,scheitelCol=0;
+                max(g,Aurora::Column,scheitelRow,scheitelCol);
+                long indexRow=0, indexCol=0;
+                //不知道对不对？？？
+                auto scheitel = scheitelRow+scheitelCol;
+                max(mean_energySum_f.block(0,0,scheitel),Aurora::Column,indexRow,indexCol);
+                size_t index = indexRow+indexCol;
+                double median_mean_energySum = Aurora::median(mean_energySum).getScalar();
+                bool flag = false; 
+                if (mean_energySum[index] > median_mean_energySum)
+                {
+                    for (size_t i = index; i > 0 ; i--)
+                    {
+                        if(mean_energySum[i]<median_mean_energySum){
+                            index = i;
+                            flag = true;
+                            break;
+                        }
+                    }
+                    if (!flag) index = 0;
+                }
+                auto hIndex_s =mean_energySum.block(0, scheitelCol, mean_energySum.getDataSize()-1)<(0.1*mean_energySum[scheitel]);
+                std::vector<int> hIndex;
+                for (size_t i = 0; i < hIndex_s.getDataSize(); i++)
+                {
+                    if (hIndex_s[i]>0)hIndex.push_back(i);
+                }
+                if (hIndex.empty()){
+                    hIndex.push_back(mean_energySum.getDataSize()-scheitel-1);
+                }
+                int index2 = hIndex[0]+scheitel-1+10;
+                //???存在没有index？？？matlab中有这一分支
+                auto sig_begin = index + start[idx] ;
+                //在这里+1 调整会matlab的index便于后续计算的正确性
+                auto sig_end = index2 + start[idx] +1;
+
+                winLength = sig_end+winLength+200-sig_begin+1;
+                auto win = generateGaussWindow(winLength);
+                double* begin = _blockedAScans.getData() + (size_t)sig_begin-1;
+                int length = winLength;
+                auto winReverse = Aurora::zeros(1,win.getDataSize());
+                std::reverse_copy(win.getData(), win.getData()+win.getDataSize(), winReverse.getData());
+                vdMulI(length, begin, 1, winReverse.getData(), 1, begin, 1);
+                winLength = sig_end;
+                win = generateGaussWindow(winLength);
+                begin = _blockedAScans.getData() ;
+                length = sig_end;
+                vdMulI(length, begin, 1, win.getData(), 1, begin, 1);
+            }
+        }
+        _blockedAScans = fft(_blockedAScans);
+        if(reflectParams::useOptPulse==1){
+            auto n = nSamples;
+            auto nHalf = round((double)n/2);
+            _blockedAScans(0,Aurora::$) = std::complex<double>{0,0};
+            _blockedAScans.setBlockComplexValue(0,nHalf,_blockedAScans.getDimSize(0)-1, std::complex<double>{0,0});
+            _blockedAScans.setBlock(0,1,nHalf-1,_blockedAScans.block(0,1,nHalf-1)*2);
+    
+    
+            _blockedAScans = ifft(_blockedAScans, n);
+    
+            _blockedAScans = abs(_blockedAScans.block(0,0,nSamples-1));
+            auto help = _blockedAScans.deepCopy();
+            help(0,Aurora::$) = 0;
+            help.setBlock(0, 1, help.getDimSize(0)-1, _blockOp(_blockedAScans, 1));
+            {
+                auto help_bbegin = help.block(0, 0, help.getDimSize(0)-2);
+                auto help_bend = help.block(0, 1, help.getDimSize(0)-1);
+                help_bend = (help_bend>0)*(help_bbegin<0);
+                auto tempBlock = Aurora::zeros(_blockedAScans.getDimSize(0),_blockedAScans.getDimSize(1),_blockedAScans.getDimSize(2));
+                tempBlock.setBlock(0, 0, tempBlock.getDimSize(0)-2, help_bend);
+                _blockedAScans = _blockedAScans*tempBlock;
+            }
+            help = _blockedAScans.deepCopy();
+            Aurora::compareSet(help,0,0,Aurora::LT);
+     
+        
+            if(reflectParams::limitNumPulsesTo>0)
+            {
+                
+                if(size(help,1)>reflectParams::limitNumPulsesTo)
+                {
+                    auto help2 = Aurora::sort(help);
+                    auto temp = repmat(help2(help2.getDimSize(0)-reflectParams::limitNumPulsesTo,Aurora::$).toMatrix(),help.getDimSize(0),1);
+                    Aurora::compareSet(help,temp,0,Aurora::LT);
+                }
+                    
+            }
+            
+            if(reflectParams::normalizePeaks)
+            {
+                Aurora::compareSet(help,0,1,Aurora::GT);
+            }
+    
+    
+            
+            help= real(ifft(fft(help)*(preComputes.sincPeak_ft.block(1,0,nAScans-1))));
+            Aurora::compareSet(_blockedAScans,0,help,Aurora::NL);
+
+        }
+        else{
+             _blockedAScans = real(ifft(_blockedAScans));
+        }
+        result.push_back(_blockedAScans);
+        
+        result.push_back(valid);
+        ascanMapValue[0] = (float)ascanMapValue[0];
+        result.push_back(ascanMapValue);
+        return result;
+    }
+
+
+    preprocessAScanRResult preprocessAScanBlockForReflection(
+    Aurora::Matrix &AscanBlock, const Aurora::Matrix &blockedMP,
+    const Aurora::Matrix &blockedSL, const Aurora::Matrix &blockedSN,
+    const Aurora::Matrix &blockedRL, const Aurora::Matrix &blockedRN,
+    const Aurora::Matrix &blockedSenderPosition,
+    const Aurora::Matrix &blockedReceiverPosition,
+    const Aurora::Matrix &blockedGain, const Aurora::Matrix &blockedChannels,
+    ExpInfo &info, const PreComputes &preComputes) 
+    {
+        preprocessAScanRResult result;
+        auto valid = Aurora::zeros(1, blockedMP.getDataSize()) ;
+        auto ascanMapValue = Aurora::zeros(1, blockedMP.getDataSize()) ;
+        // auto dOutMaxVal = Aurora::zeros(1, blockedMP.getDataSize()) ;
+        // auto dOutStdVal = Aurora::zeros(1, blockedMP.getDataSize()) ;
+        // auto dOutMeanVal = Aurora::zeros(1, blockedMP.getDataSize()) ;
+        // auto dOutSnr = Aurora::zeros(1, blockedMP.getDataSize()) ;
+        // auto dOutSnrPass = Aurora::zeros(1, blockedMP.getDataSize()) ;
+
+        switch (reflectParams::version) {
+            case 1:
+                std::cerr<<"error USCT II only!!"<<std::endl;
+                break;
+            case 2:{
+                #pragma omp parallel for num_threads(32)
+                for (size_t i = 0; i < blockedMP.getDataSize(); i++)
+                {
+                    auto signalPResult = performSignalProcessing(
+                        AscanBlock(Aurora::$, i).toMatrix(),
+                        blockedSL(Aurora::$, i).toMatrix(),
+                        blockedRL(Aurora::$, i).toMatrix(),
+                        blockedSenderPosition(Aurora::$, i).toMatrix(),
+                        blockedReceiverPosition(Aurora::$, i).toMatrix(),
+                        blockedGain(Aurora::$, i).toMatrix(),
+                        blockedChannels(Aurora::$, i).toMatrix(), info,
+                        preComputes);
+                    AscanBlock(Aurora::$,i)    = signalPResult[0];
+                    valid(Aurora::$,i)         = signalPResult[1];
+                    ascanMapValue(Aurora::$,i) = signalPResult[2];
+
+                // dOutMaxVal[idx] = dOut.maxVal;
+                // dOutStdVal[idx] = dOut.stdVal;
+                // dOutMeanVal[idx] = dOut.meanVal;
+                // dOutSnr[idx] = dOut.snr;
+                // dOutSnrPass[idx] = dOut.snrPass;
+                }
+            }
+        }
+        if (reflectParams::findDefects==1)
+        {
+            auto highNoiseScore = Aurora::zeros(1, ascanMapValue.getDataSize());
+            size_t count = 0;
+            for (size_t i = 0; i < valid.getDataSize(); i++)
+            {
+                if (valid[i])highNoiseScore[count++] = ascanMapValue[i];
+            }
+            highNoiseScore.forceReshape(1, count, 1);
+            double highNoiseScoreMean = mean(highNoiseScore).getScalar();
+            double highNoiseScoreStd = Aurora::std(highNoiseScore).getScalar();
+            double treshold99= (mean(highNoiseScore)+2.596*Aurora::std(highNoiseScore)).getScalar(); 
+            Aurora::compareSet(valid,ascanMapValue,treshold99,0,Aurora::GT);
+        // debugOutput.blockStatisticValue = ascanMapValue;
+        // debugOutput.blockStatisticUsedData = ascanMapValue > treshold99;
+        // debugOutput.maxVal = dOutMaxVal;
+        // debugOutput.stdVal = dOutStdVal;
+        // debugOutput.meanVal = dOutMeanVal;
+        // debugOutput.snr = dOutSnr;
+        // debugOutput.snrPass = dOutSnrPass;
+        }
+        result.AscanBlock = AscanBlock;
+        result.usedData = valid;
+        return result;
+    }
+} // namespace Recon
diff --git a/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.h b/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.h
new file mode 100644
index 0000000..d3ae9ce
--- /dev/null
+++ b/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.h
@@ -0,0 +1,26 @@
+#ifndef __PREPROCESSASCANBLOCKFORREFLECTION_H__
+#define __PREPROCESSASCANBLOCKFORREFLECTION_H__
+#include "Matrix.h"
+#include "common/getMeasurementMetaData.h"
+#include <cstddef>
+namespace Recon {
+    struct preprocessAScanRResult {
+    Aurora::Matrix AscanBlock;
+    Aurora::Matrix usedData;
+    };
+    struct ExpInfo{
+        size_t Wavelength;
+        size_t expectedAScanLength;
+    };
+
+    preprocessAScanRResult preprocessAScanBlockForReflection(
+        Aurora::Matrix &AscanBlock, const Aurora::Matrix &mpBlock,
+        const Aurora::Matrix &slBlock, const Aurora::Matrix &snBlock,
+        const Aurora::Matrix &rlBlock, const Aurora::Matrix &rnBlock,
+        const Aurora::Matrix &senderPositionBlock,
+        const Aurora::Matrix &receiverPositionBlock,
+        const Aurora::Matrix &gainBlock, const Aurora::Matrix &channelBlock,
+        ExpInfo& info, const PreComputes &preComputes);
+}
+
+#endif // __PREPROCESSASCANBLOCKFORREFLECTION_H__
\ No newline at end of file
diff --git a/test/Reconstruction_Test.cpp b/test/Reconstruction_Test.cpp
index de2a39b..8e1ccd2 100644
--- a/test/Reconstruction_Test.cpp
+++ b/test/Reconstruction_Test.cpp
@@ -8,6 +8,7 @@
 #include "Matrix.h"
 #include "Sparse.h"
 #include "config/config.h"
+#include "reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.h"
 #include "transmissionReconstruction/reconstruction/buildMatrix/DGradient.h"
 #include "transmissionReconstruction/reconstruction/buildMatrix/FMM.h"
 #include "transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h"
@@ -46,16 +47,48 @@ protected:
 };
 
 TEST_F(Reconstruction_Test, determineOptimalPulse) {
-    Recon::reflectParams::imageResolution = 8.6381e-04;
-    Recon::reflectParams::optPulseFactor = 24;
-    auto result = Recon::determineOptimalPulse(1.0000e-07,3000);
-    EXPECT_EQ(3000, result.getDataSize());
-    ASSERT_DOUBLE_AE(result[2], 0.0025);
-    ASSERT_DOUBLE_AE(result[4], 0.0078);
+    Recon::reflectParams::imageResolution = 8.925316499483377e-04;
+    Recon::reflectParams::optPulseFactor = 48;
+    auto result = Recon::determineOptimalPulse(1.0000e-07,4096);
+    EXPECT_EQ(4096, result.getDataSize());
+    MatlabReader m2("/home/krad/TestData/sincPeakft.mat");
+    auto f1 = m2.read("sincPeak_ft");
+    for(size_t i=0; i<f1.getDataSize(); ++i)
+    {
+        EXPECT_DOUBLE_AE(f1[i], 0)<<"index:"<<i;
+    }
+}
 
-    // for(size_t i=0; i<result.getDataSize(); ++i)
+TEST_F(Reconstruction_Test, preprocessAScanBlockForReflection) {
+    Recon::initalizeConfig();
+    MatlabReader m2("/home/krad/TestData/preprocessRefC.mat");
+    auto blockedAScans = m2.read("blockedAScans");
+    auto blockedMP = m2.read("blockedMP");
+    auto blockedSL = m2.read("blockedSL");
+    auto blockedSN = m2.read("blockedSN");
+    auto blockedRL = m2.read("blockedRL");
+    auto blockedRN = m2.read("blockedRN");
+    auto senderPositionBlock = m2.read("blockedSenderPosition");
+    auto receiverPositionBlock = m2.read("blockedReceiverPosition");
+    auto channelBlock = m2.read("blockedChannels");
+    auto gainBlock = m2.read("blockedGain");
+    Recon::ExpInfo info;
+    info.expectedAScanLength= 4096;
+    info.Wavelength=3;
+    Recon::PreComputes preComputes;
+    preComputes.measuredCEUsed = true;
+    preComputes.offset = 5.2000e-07;
+    preComputes.timeInterval = 1.0000e-07;
+    preComputes.matchedFilter =  m2.read("matchedFilter");
+    MatlabReader m("/home/krad/TestData/sincPeakft.mat");
+    preComputes.sincPeak_ft = m.read("sincPeak_ft");
+    auto result = Recon::preprocessAScanBlockForReflection(blockedAScans, blockedMP, blockedSL, blockedSN, blockedRL,
+     blockedRN,senderPositionBlock, receiverPositionBlock, gainBlock, channelBlock, info, preComputes);
+
+    size_t size = result.AscanBlock.getDataSize();
+    // for(size_t i=0; i<result.AscanBlock.getDataSize(); ++i)
     // {
-    //     ASSERT_DOUBLE_AE(f1[i], result.outSOS[i])<<"index:"<<i;
+    //     ASSERT_DOUBLE_AE(f1[i], result[i])<<"index:"<<i;
     // }
 }