From 5e003c734823a792753a0cd2f23753a4b4fd4965 Mon Sep 17 00:00:00 2001
From: kradchen <chenjuijun@realcan.cn>
Date: Fri, 9 Jun 2023 17:55:40 +0800
Subject: [PATCH] Add determineOptimalPulse function

---
 .../preprocessData/determineOptimalPulse.cpp  | 67 +++++++++++++++++++
 .../preprocessData/determineOptimalPulse.h    | 11 +++
 test/Reconstruction_Test.cpp                  | 43 ++++++++++++
 3 files changed, 121 insertions(+)
 create mode 100644 src/reflectionReconstruction/preprocessData/determineOptimalPulse.cpp
 create mode 100644 src/reflectionReconstruction/preprocessData/determineOptimalPulse.h

diff --git a/src/reflectionReconstruction/preprocessData/determineOptimalPulse.cpp b/src/reflectionReconstruction/preprocessData/determineOptimalPulse.cpp
new file mode 100644
index 0000000..65ef3c3
--- /dev/null
+++ b/src/reflectionReconstruction/preprocessData/determineOptimalPulse.cpp
@@ -0,0 +1,67 @@
+#include "determineOptimalPulse.h"
+#include "Function1D.h"
+#include "Function2D.h"
+#include "Function3D.h"
+#include "config/config.h"
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <cstring>
+#include <omp.h>
+namespace Recon {
+
+    Aurora::Matrix determineOptimalPulse(double timeInterval, size_t expectedAScanLength)
+    {
+        double imgResolution = reflectParams::imageResolution;
+        int optPulseFactor = reflectParams::optPulseFactor;
+        auto minOptPulse = ceil(8*imgResolution*(1/timeInterval)/1500);
+
+        if(optPulseFactor==-1)
+        {
+            optPulseFactor = minOptPulse;
+            reflectParams::optPulseFactor = -1;
+            printf("Optimal pulse automatically set to %f \n",minOptPulse);
+        }
+        else if(optPulseFactor< minOptPulse){
+            printf("WARNING: optimal pulse too small for resolution");
+        }
+
+        auto desSamplingFreq = (1/timeInterval)*0.5;
+        double sincFact = optPulseFactor;
+        double sincLength = round(11*sincFact/16);
+        sincLength = 3*sincLength; 
+        double begin =-sincLength/2+(timeInterval)/2;
+        double end = sincLength/2;
+        int length = end - begin + 1;
+        auto sincT = Aurora::zeros(1,length);
+        for (int j = 0; begin <= end; begin++, j++)
+        {
+            sincT[j] = begin*timeInterval;
+        }
+
+        sincFact = sincFact/2;
+        auto sincPeak = 2*sqrt(M_PI)*std::pow((desSamplingFreq/sincFact),3)*(1-2*((M_PI*desSamplingFreq/sincFact*sincT)^2))*exp(-((M_PI*desSamplingFreq/sincFact*sincT)^2));
+
+        sincPeak = sincPeak/Aurora::max(sincPeak);
+        auto sincPeak_len=sincPeak.getDataSize();
+
+        Aurora::padding(sincPeak, expectedAScanLength-1, 0);
+        // sincPeak = Aurora::transpose(sincPeak);
+        size_t offset = floor((double)sincPeak_len/2);
+        //cicshift
+        // #pragma omp  parallel for
+        for (size_t i = 0; i < sincPeak.getDimSize(1); i++)
+        {
+            double *temp = new double[offset]{0};
+            double * beginPtr = sincPeak.getData()+i*sincPeak.getDimSize(0);
+            double * endPtr = beginPtr + offset;
+            double * col_end = beginPtr + (i+1)*sincPeak.getDimSize(0);
+            std::copy(beginPtr,endPtr,temp);
+            std::copy(endPtr,col_end,beginPtr);
+            std::copy(temp,temp+offset,col_end-offset);
+        }
+        auto sincPeak_ft=fft(sincPeak);
+        return sincPeak_ft;
+    }
+}
+
diff --git a/src/reflectionReconstruction/preprocessData/determineOptimalPulse.h b/src/reflectionReconstruction/preprocessData/determineOptimalPulse.h
new file mode 100644
index 0000000..d4a695f
--- /dev/null
+++ b/src/reflectionReconstruction/preprocessData/determineOptimalPulse.h
@@ -0,0 +1,11 @@
+#ifndef __DETERMINEOPTIMALPULSE_H__
+#define __DETERMINEOPTIMALPULSE_H__
+
+#include "Matrix.h"
+namespace Recon {
+Aurora::Matrix determineOptimalPulse(double timeInterval,
+                                     size_t expectedAScanLength);
+
+}
+
+#endif // __DETERMINEOPTIMALPULSE_H__
\ No newline at end of file
diff --git a/test/Reconstruction_Test.cpp b/test/Reconstruction_Test.cpp
index ece1767..de2a39b 100644
--- a/test/Reconstruction_Test.cpp
+++ b/test/Reconstruction_Test.cpp
@@ -4,6 +4,7 @@
 
 #include "Function1D.h"
 #include "MatlabReader.h"
+#include "MatlabWriter.h"
 #include "Matrix.h"
 #include "Sparse.h"
 #include "config/config.h"
@@ -13,6 +14,8 @@
 #include "transmissionReconstruction/reconstruction/reconstruction.h"
 #include "transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h"
 
+#include "reflectionReconstruction/preprocessData/determineOptimalPulse.h"
+
 
 
 
@@ -42,6 +45,20 @@ 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);
+
+    // for(size_t i=0; i<result.getDataSize(); ++i)
+    // {
+    //     ASSERT_DOUBLE_AE(f1[i], result.outSOS[i])<<"index:"<<i;
+    // }
+}
+
 TEST_F(Reconstruction_Test, reconstructArt) {
     MatlabReader m("/home/sun/testData/reconstructArt.mat");
     auto data  = m.read("data");
@@ -261,3 +278,29 @@ TEST_F(Reconstruction_Test,traceStraightRayBresenham){
     EXPECT_DOUBLE_AE(result.path[7],2);
     EXPECT_EQ(6, result.pathLen);
 }
+
+TEST_F(Reconstruction_Test,callTval3){
+    TVALOptions opt;
+    opt.nonneg = false;
+    opt.bent = false;
+    opt.tol = 1E-10;
+    opt.maxit = 50;
+    opt.TVnorm = 2;
+    opt.disp = false;
+    opt.mu0 = 100;
+    opt.mu = 100;
+    opt.beta = 1;
+    opt.beta0 = 1;
+    MatlabReader m("/home/sun/testData/buildMatrixTest.mat");
+    auto j1 = m.read("j1");
+    auto i1 = m.read("i1");
+    auto s1 = m.read("s1");
+    Aurora::Sparse M(i1-1,j1-1,s1,734989,1196032);
+    MatlabReader m2("/home/sun/testData/tval3gpu3d.mat");
+    auto b = m2.read("b");
+    auto dims = Aurora::Matrix::fromRawData(new double[3]{128,128,73}, 1, 3);
+    auto result = Recon::callTval3(M, b, dims, 0,opt);
+    auto outSOS = Recon::slownessToSOS(result, 1.498206569328594e+03) ;
+    MatlabWriter w2("/home/krad/transmissionSOS111.mat");
+    w2.write(outSOS, "SOS");
+}