feat: Add logic block log to performSignalProcessing in transmission reconstruction

feat: Add logic block log to performSignalProcessing in reflect reconstruction
feat: A CPU Performance Test tool branch, which remove cuda & outpout logic.
2024-11-06 13:57:11 +08:00 · 2024-11-06 11:36:52 +08:00 · 2024-11-04 10:55:21 +08:00 · 2024-10-22 10:25:48 +08:00
26 changed files with 317 additions and 1161 deletions
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -4,80 +4,23 @@ project(UR)

 set(CMAKE_CXX_STANDARD 14)
 set(CMAKE_INCLUDE_CURRENT_DIR ON)
-add_definitions(-DUSE_CUDA)
-set(CMAKE_CUDA_COMPILER /usr/local/cuda/bin/nvcc)
-enable_language(CUDA)
-find_package(CUDAToolkit REQUIRED)
-set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -g")

+set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -g")
+set(MKL_LINK static)
 find_package(Aurora REQUIRED)
 find_package(Parser REQUIRED)
-find_package(URDepends REQUIRED)
-find_package(DCMTK REQUIRED)
-set(DCMTK_WITH_JPEG OFF)
 file(GLOB_RECURSE  cpp_files ./src/*.cpp)
-file(GLOB_RECURSE  cu_files ./src/*.cu)
 file(GLOB_RECURSE  cxx_files ./src/*.cxx)
 file(GLOB_RECURSE  header_files ./src/*.h)
-add_executable(UR  ${cpp_files} ${cu_files} ${cxx_files} ${header_files} ${Aurora_Source} ${Aurora_Source_Cu} ./src/Aurora.cu ${Aurora_DIR}/src/CudaMatrixPrivate.cu)
+add_executable(UR  ${cpp_files} ${cxx_files} ${header_files} ${Aurora_Source})

 #target_compile_options(UR PUBLIC ${Aurora_Complie_Options} "-march=native")
 target_include_directories(UR PUBLIC ./src/)
 target_include_directories(UR PUBLIC ${Aurora_INCLUDE_DIRS})
-target_include_directories(UR PUBLIC ${DCMTK_INCLUDE_DIRS})
 target_include_directories(UR PUBLIC ${Parser_INCLUDE_DIRS})
-target_include_directories(UR PUBLIC ${URDepends_INCLUDES_DIRS})
-target_include_directories(UR PUBLIC ${DCMTK_INCLUDE_DIRS})
+
 target_link_libraries(UR PUBLIC ${Aurora_Libraries})
-target_link_libraries(UR PUBLIC dcmdata)
+
 target_link_libraries(UR PUBLIC matio)
-target_link_libraries(UR PUBLIC Parser)
-
-target_link_libraries(UR PUBLIC URDepends::TransDetection)
-target_link_libraries(UR PUBLIC URDepends::eikonal)
-target_link_libraries(UR PUBLIC URDepends::TVALGPU)
-target_link_libraries(UR PUBLIC URDepends::SaftTofi)
-target_link_libraries(UR PUBLIC dcmdata)
-
-
-target_include_directories(UR PRIVATE ./src /usr/local/cuda/include) 
-set_target_properties(UR PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
-#target_compile_options(UR PRIVATE ${Aurora_Complie_Options} "-march=native")
-# target_compile_options(UR PRIVATE $<$<COMPILE_LANGUAGE:CUDA>:
-#                        --use_fast_math
-#                        --ptxas-options=-v
-#                        -arch sm_75
-#                        >)
-target_compile_options(UR PRIVATE $<$<COMPILE_LANGUAGE:CUDA>:
-    -arch=sm_75 --expt-extended-lambda 
->)
-
-target_link_libraries(UR PRIVATE ${CUDA_RUNTIME_LIBRARY} CUDA::cufft CUDA::cudart)
-
-# target_link_libraries(UR PUBLIC URDepends::SaftATT)
-find_package(GTest  REQUIRED)
-INCLUDE_DIRECTORIES(${GTEST_INCLUDE_DIRS})
-
-enable_testing()
-file(GLOB_RECURSE  test_cpp test/*.cpp)
-add_executable(UR_Test ${cpp_files} ${header_files} ${Aurora_Source} ${test_cpp} )
-target_include_directories(UR_Test PUBLIC ./test/ ./src/)
-# target_compile_options(UR_Test PUBLIC ${Aurora_Complie_Options} "-march=native")
-target_include_directories(UR_Test PUBLIC ${DCMTK_INCLUDE_DIRS})
-target_include_directories(UR_Test PUBLIC ${Aurora_INCLUDE_DIRS})
-target_include_directories(UR_Test PUBLIC ${Parser_INCLUDE_DIRS})
-target_include_directories(UR_Test PUBLIC ${URDepends_INCLUDES_DIRS})
-# 必须写前面不然容易出问题
-target_link_libraries(UR_Test PUBLIC ${GTEST_BOTH_LIBRARIES})
-target_link_libraries(UR_Test PUBLIC ${Aurora_Libraries})
-target_link_libraries(UR_Test PUBLIC dcmdata)
-target_link_libraries(UR_Test PUBLIC matio)
-# target_link_libraries(UR_Test PUBLIC ${Parser_Libraries})
-target_link_libraries(UR_Test PUBLIC Parser)
-target_link_libraries(UR_Test PUBLIC URDepends::TransDetection)
-target_link_libraries(UR_Test PUBLIC URDepends::eikonal)
-target_link_libraries(UR_Test PUBLIC URDepends::TVALGPU)
-target_link_libraries(UR_Test PUBLIC URDepends::SaftTofi)
-# target_link_libraries(UR_Test PUBLIC URDepends::SaftATT)
-
-gtest_discover_tests(UR_Test)
+target_link_libraries(UR PUBLIC pthread)
+target_link_libraries(UR PUBLIC Parser)
--- a/src/common/DICOMExporter.cpp
+++ b/src/common/DICOMExporter.cpp
@@ -1,260 +0,0 @@
-#include "DICOMExporter.h"
-
-#include "fileHelper.h"
-#include "config/config.h"
-#include "log/log.h"
-
-#include <cstddef>
-#include <cstdio>
-#include <dcmtk/dcmdata/dcdeftag.h>
-#include <dcmtk/dcmdata/dcfilefo.h>
-#include <dcmtk/dcmdata/dcuid.h>
-#include <dcmtk/ofstd/offile.h>
-#include <string>
-#include <sys/types.h>
-#include "dcmtk/dcmdata/dctk.h"
-
-using namespace std;
-namespace  {
-    float getSpacing(float * startPoint, float* endPoint, int* imageXYZ, int index){
-    return (endPoint[index]-startPoint[index])/(float)(imageXYZ[index])*1000.0;
-}
-
-    void initialDicomFile(DcmDataset* dataset,DcmMetaInfo* metaInfo)
-    {
-        metaInfo->putAndInsertString(DCM_MediaStorageSOPClassUID,
-                                        UID_CTImageStorage);
-        dataset->putAndInsertString(DCM_SOPClassUID, UID_CTImageStorage);
-
-        dataset->putAndInsertString(DCM_Modality, "CT");
-        dataset->putAndInsertString(DCM_TransferSyntaxUID,
-                                        UID_LittleEndianImplicitTransferSyntax);
-        dataset->putAndInsertString(DCM_ImageType, "ORIGINAL\\PRIMARY");
-
-        dataset->putAndInsertString(DCM_Manufacturer, "EQ9");
-        dataset->putAndInsertString(DCM_PatientOrientation, "L\\H");
-        dataset->putAndInsertString(DCM_PatientPosition, "HFP");
-
-        dataset->putAndInsertString(DCM_ImageOrientationPatient, "1\\0\\0\\0\\0\\1");
-
-        dataset->putAndInsertString(DCM_StudyDescription, "USCT BREAST");
-        dataset->putAndInsertString(DCM_ManufacturerModelVersion, "USCT3DV3");
-
-        dataset->putAndInsertString(DCM_SamplesPerPixel, "1");
-        dataset->putAndInsertString(DCM_PixelRepresentation, "0");
-        dataset->putAndInsertString(DCM_BitsAllocated, "16");
-        dataset->putAndInsertString(DCM_BitsStored, "16");
-        dataset->putAndInsertString(DCM_HighBit, "15");
-        dataset->putAndInsertString(DCM_PhotometricInterpretation, "MONOCHROME2");
-
-        dataset->putAndInsertString(DCM_KVP, "0");
-
-    }
-    void initialDate(DcmDataset* dataset,const std::string& mStudyDate )
-    {
-        dataset->putAndInsertString(DCM_StudyDate, mStudyDate.data());
-        OFDateTime currentDateTime;
-        currentDateTime.setCurrentDateTime();
-        OFString date;
-        currentDateTime.getDate().getISOFormattedDate(date,false);
-        dataset->putAndInsertString(DCM_SeriesDate, date.data());
-        dataset->putAndInsertString(DCM_AcquisitionDate, date.data());
-        dataset->putAndInsertString(DCM_ContentDate, date.data());
-    }
-
-    void initialSeriesTime(DcmDataset* dataset){
-                // 获取当前日期和时间（UTC时间）
-        OFDateTime currentDateTime;
-        currentDateTime.setCurrentDateTime();
-        OFString time;
-        currentDateTime.getTime().getISOFormattedTime(time,true,false,false,false);
-        dataset->putAndInsertString(DCM_SeriesTime, time.data());
-        dataset->putAndInsertString(DCM_AcquisitionTime, time.data());
-        dataset->putAndInsertString(DCM_ContentTime, time.data());
-    }
-    const char* SeriesDescriptions[3]={"REFL","SOS","ATT"};
-    const char* getSeriesDescription(int type){
-        return SeriesDescriptions[type];
-    }
-}
-
-
-namespace Recon
-{
-    DICOMExporter::DICOMExporter(const PatientData& aPatientData)
-    {
-        mPatientData = aPatientData;
-        //prepare Study
-        dcmGenerateUniqueIdentifier(mStudyInstanceUID,SITE_STUDY_UID_ROOT);
-
-        OFDateTime currentDateTime;
-        currentDateTime.setCurrentDateTime();
-        OFString date;
-        currentDateTime.getDate().getISOFormattedDate(date,false);
-        OFString time;
-        currentDateTime.getTime().getISOFormattedTime(time,true,false,false,false);
-        mStudyDate = date.data();
-        mStudyTime = time.data();
-    }
-    
-    DICOMExporter::DICOMExporter()
-    {
-                //prepare Study
-        dcmGenerateUniqueIdentifier(mStudyInstanceUID,SITE_STUDY_UID_ROOT);
-
-        OFDateTime currentDateTime;
-        currentDateTime.setCurrentDateTime();
-        OFString date;
-        currentDateTime.getDate().getISOFormattedDate(date,false);
-        OFString time;
-        currentDateTime.getTime().getISOFormattedTime(time,true,false,false,false);
-        mStudyDate = date.data();
-        mStudyTime = time.data();
-    }
-
-    DICOMExporter::~DICOMExporter()
-    {
-    }
-    void DICOMExporter::setExportBasePath(const std::string & path)
-    {
-        mBasePath = path;
-    }
-    
-    void DICOMExporter::exportDICOM(Aurora::Matrix aMatrix, ImageType type){
-        int XYZ [3] = {aMatrix.getDimSize(0),aMatrix.getDimSize(1), aMatrix.getDimSize(2)};
-        string path = mBasePath+"/"+::getSeriesDescription(type)+"/";
-        if (!exists(path)){
-            if(!mkdir(path)){
-                RECON_ERROR("Fail to create dir with path {0}\r\n",path);
-            }
-        }
-        exportDCM(path, reflectParams::imageStartpoint.getData(),
-            reflectParams::imageEndpoint.getData(), XYZ, aMatrix.getData(), type);
-        RECON_INFO("Save DICOM to Path:{0}, type:{1}==================>",path,(int)type);
-    }
-
-    void DICOMExporter::exportDCM(string path, float * startPoint, float* endPoint, int* imageXYZ,float* data, int type)
-    {
-        long Rows = imageXYZ[1], Cols = imageXYZ[0] ,Slices = imageXYZ[2];
-        DcmFileFormat dcmFile;
-        DcmDataset* dataset = dcmFile.getDataset();
-        DcmMetaInfo* metaInfo = dcmFile.getMetaInfo();
-
-        initialDicomFile(dataset,metaInfo);
-
-        dataset->putAndInsertString(DCM_StudyInstanceUID, mStudyInstanceUID);
-
-        dcmGenerateUniqueIdentifier(mSeriesInstanceUID,SITE_SERIES_UID_ROOT);
-        dataset->putAndInsertString(DCM_SeriesInstanceUID, mSeriesInstanceUID);
-        dataset->putAndInsertString(DCM_FrameOfReferenceUID, mSeriesInstanceUID);
-
-        ::initialDate(dataset, mStudyDate);
-
-        dataset->putAndInsertString(DCM_StudyTime, mStudyTime.data());
-
-        ::initialSeriesTime(dataset);
-
-        dataset->putAndInsertString(DCM_InstitutionName, mPatientData.getInstituationName().empty()?
-            "InstName":mPatientData.getInstituationName().data());
-        dataset->putAndInsertString(DCM_InstitutionAddress, mPatientData.getInstituationAddress().empty()?
-            "default addr":mPatientData.getInstituationAddress().data());
-        dataset->putAndInsertString(DCM_ReferringPhysicianName,mPatientData.getReferringPhysicianName().empty()?
-            "ReferringPhysician": mPatientData.getOperatorName().data());
-        dataset->putAndInsertString(DCM_PatientName,  mPatientData.getPatientName().empty()?
-            "TestPatient":mPatientData.getPatientName().data());
-        dataset->putAndInsertString(DCM_PatientSex, mPatientData.getPatientSex().data());
-        dataset->putAndInsertString(DCM_PatientBirthDate, mPatientData.getPatientBirthDate().data());
-        dataset->putAndInsertString(DCM_PatientID, mPatientData.getPatientID().empty()?
-            "TestPatID":mPatientData.getPatientID().data());
-        dataset->putAndInsertString(DCM_Laterality, mPatientData.getLaterality().data());
-
-        dataset->putAndInsertString(DCM_StudyID, "0");
-        dataset->putAndInsertString(DCM_SeriesNumber, to_string(type).data());
-        dataset->putAndInsertString(DCM_SeriesDescription, ::getSeriesDescription(type));
-        // dataset->putAndInsertString(DCM_ProtocolName, "?");
-
-        // initial spacing and position
-        float spacing[3]{.0,.0,.0};
-        spacing[0] = getSpacing(startPoint,endPoint,imageXYZ,0);
-        spacing[1] = getSpacing(startPoint,endPoint,imageXYZ,1);
-        spacing[2] = getSpacing(startPoint,endPoint,imageXYZ,2);
-
-        float originPosition[3] = {
-            endPoint[1]*(float)1000.0,
-            endPoint[2]*(float)1000.0,
-            endPoint[0]*(float)1000.0,
-        };
-        float originLocation =-endPoint[2]*1000.0;
-        dataset->putAndInsertString(DCM_SliceThickness, to_string(spacing[2]).data());
-        dataset->putAndInsertUint16(DCM_Rows, Rows);
-        dataset->putAndInsertUint16(DCM_Columns, Cols);
-        string spacingsBP = to_string(spacing[0])+"\\"+to_string(spacing[1]);
-        dataset->putAndInsertString(DCM_PixelSpacing, spacingsBP.data());
-
-        // initial data, wwwl , Slope,Intercept 
-        float Slope = 1, Intercept = 0 ;
-        size_t size = Rows*Cols*Slices;
-        ushort* udata = new ushort[size]{0};
-        float min = data[0];
-        float max = data[0];
-        for (size_t i = 0; i < size; i++)
-        {
-            if (type == 0) data[i] = 1000*data[i];
-            if (min > data[i]) min = data[i];
-            if (max < data[i]) max = data[i];
-        }
-
-        long windowCenter = min + (max-min)/2, windowWidth = max-min;
-        float dSlope = 1.0/(pow(2,16)/windowWidth);
-        float dIntercept  = min;
-        for (size_t i = 0; i < size; i++)
-        {
-            udata[i] = (ushort)((data[i] - min)/dSlope);
-        }
-        if (type !=0){
-            Slope = dSlope;
-            Intercept = dIntercept;
-        }
-        else{
-            windowCenter = (windowCenter-min)/dSlope;
-            windowWidth = (windowWidth-min)/dSlope;
-        }
-        dataset->putAndInsertString(DCM_WindowCenter, to_string(windowCenter).data());
-        dataset->putAndInsertString(DCM_WindowWidth, to_string(windowWidth).data());
-        dataset->putAndInsertString(DCM_RescaleSlope, to_string(Slope).data());
-        dataset->putAndInsertString(DCM_RescaleIntercept, to_string(Intercept).data());
-
-        size_t sliceSize = Rows*Cols;
-        for (size_t i = 0; i < Slices; i++)
-        {
-            dataset->putAndInsertString(DCM_AccessionNumber, to_string(i).data());
-            dataset->putAndInsertString(DCM_InstanceNumber, to_string(i).data());
-            string pos = to_string(originPosition[1] ) + "\\" +
-                        to_string(originPosition[2]+ i * spacing[2]) + "\\" +
-                        to_string(originPosition[0]);
-            dataset->putAndInsertString(DCM_ImagePositionPatient, pos.data());
-            string loc = to_string(originLocation + i * spacing[2]);
-            dataset->putAndInsertString(DCM_SliceLocation, loc.data());
-
-            // set SOPInstanceUID
-            char SOP_UID[100] = {0};
-            dcmGenerateUniqueIdentifier(SOP_UID, SITE_INSTANCE_UID_ROOT);
-            metaInfo->putAndInsertString(DCM_MediaStorageSOPInstanceUID, SOP_UID);
-            dataset->putAndInsertString(DCM_SOPInstanceUID, SOP_UID);
-
-            dataset->putAndInsertUint16Array(DCM_PixelData, udata + sliceSize * i,
-                                            sliceSize);
-            string path2 = path+to_string(i)+".dcm";
-            OFFilename file ;
-            file = path2.data();
-            auto status =dcmFile.saveFile(file, EXS_LittleEndianExplicit);
-            if(status.bad())
-            {
-                RECON_ERROR("Save dicom file {0} fail, beacuse of {1}\r\n",path2.data(),status.text());
-            }
-            else{
-                RECON_INFO("Save dicom file {0} \r\n",path2.data());
-            }
-        }
-    }
-}
--- a/src/common/DICOMExporter.h
+++ b/src/common/DICOMExporter.h
@@ -1,41 +0,0 @@
-#ifndef __DICOMEXPORT_H__
-#define __DICOMEXPORT_H__
-#include <string>
-
-#include "Matrix.h"
-
-#include "Parser/Data/PatientData.h"
-namespace Recon {
-
-class DICOMExporter
-{
-public:
-    enum ImageType{
-        REFL,
-        SOS,
-        ATT
-    };
-    explicit DICOMExporter(const PatientData& aPatientData);
-    //仅仅用于测试！！！
-    DICOMExporter();
-    DICOMExporter(DICOMExporter &&) = default;
-    DICOMExporter(const DICOMExporter &) = default;
-    DICOMExporter &operator=(DICOMExporter &&) = default;
-    DICOMExporter &operator=(const DICOMExporter &) = default;
-    ~DICOMExporter();
-    void setExportBasePath(const std::string & path);
-    void exportDICOM(Aurora::Matrix aMatrix, ImageType type);
-private:
-    std::string mBasePath;
-    std::string mStudyDate;
-    std::string mStudyTime;
-
-    char mStudyInstanceUID[100]={0};
-    char mSeriesInstanceUID[100]={0};
-    PatientData mPatientData;
-    void exportDCM(std::string path, float * startPoint, float* endPoint, int* imageXYZ,float* data, int type);
-};
-
-
-}
-#endif // __DICOMEXPORT_H__
--- a/src/common/dataBlockCreation/getAScanBlockPreprocessed.cpp
+++ b/src/common/dataBlockCreation/getAScanBlockPreprocessed.cpp
@@ -90,70 +90,4 @@ AscanBlockPreprocessed Recon::getAscanBlockPreprocessed(Parser* aParser, const A
    }
 //printTime();    
    return result;
-}
-
-
-AscanBlockPreprocessedCuda Recon::getAscanBlockPreprocessedCuda(Parser* aParser, const Aurora::Matrix& aMp, const Aurora::Matrix& aSl, const Aurora::Matrix& aSn, 
-                                                         const Aurora::Matrix& aRl, const Aurora::Matrix& aRn, GeometryInfo& aGeom, const MeasurementInfo& aMeasInfo,
-                                                         bool aApplyFilter, bool aTransReco)
-{
-//std::cout<<"strart"<<std::endl;
-//printTime();
-//550ms
-    AscanBlockPreprocessedCuda result;
-    AscanBlock ascanBlock = getAscanBlock(aParser, aMp, aSl, aSn, aRl, aRn);
-//printTime();
-//300ms
-    result.ascanBlockPreprocessed = ascanBlock.ascanBlock.toDeviceMatrix();
-    result.gainBlock = ascanBlock.gainBlock.toDeviceMatrix();
-    result.mpBlock = ascanBlock.mpBlock;
-    result.rlBlock = ascanBlock.rlBlock;
-    result.rnBlock = ascanBlock.rnBlock;
-    result.slBlock = ascanBlock.slBlock;
-    result.snBlock = ascanBlock.snBlock;
-    GeometryBlock geometryBlock = blockingGeometryInfos(aGeom, ascanBlock.rnBlock, ascanBlock.rlBlock, ascanBlock.snBlock, ascanBlock.slBlock, ascanBlock.mpBlock);
-//printTime();
-//3ms    
-    result.receiverPositionBlock = geometryBlock.receiverPositionBlock;
-    result.senderPositionBlock = geometryBlock.senderPositionBlock;
-    if(aApplyFilter)
-    {
-        Matrix usedData;
-        if(aTransReco)
-        {
-            usedData = filterTransmissionData(ascanBlock.slBlock, ascanBlock.snBlock, ascanBlock.rlBlock, ascanBlock.rnBlock,
-                                              aGeom.sensData, geometryBlock.senderNormalBlock, geometryBlock.receiverNormalBlock);
-        }
-        else
-        {
-            usedData = filterReflectionData(geometryBlock.receiverPositionBlock, geometryBlock.senderPositionBlock, geometryBlock.senderNormalBlock, reflectParams::constrictReflectionAngles);
-        }
-//printTime();
-//40ms
-        CudaMatrix usedDataDevice = usedData.toDeviceMatrix();
-        result.ascanBlockPreprocessed = valid(result.ascanBlockPreprocessed, usedDataDevice);
-        result.mpBlock = removeDataFromArrays(ascanBlock.mpBlock, usedData);
-        result.slBlock = removeDataFromArrays(ascanBlock.slBlock, usedData);
-        result.snBlock = removeDataFromArrays(ascanBlock.snBlock, usedData);
-        result.rlBlock = removeDataFromArrays(ascanBlock.rlBlock, usedData);
-        result.rnBlock = removeDataFromArrays(ascanBlock.rnBlock, usedData);
-
-        result.senderPositionBlock = removeDataFromArrays(geometryBlock.senderPositionBlock, usedData);
-        result.receiverPositionBlock = removeDataFromArrays(geometryBlock.receiverPositionBlock, usedData);
-        result.gainBlock = valid(result.gainBlock, usedDataDevice);
-//printTime();
-//10ms
-    }
-
-    if (ascanBlock.ascanBlock.getDataSize() > 0)
-    {
-        result.ascanBlockPreprocessed = preprocessAscanBlockCuda(result.ascanBlockPreprocessed, aMeasInfo);
-    }
-    // else
-    // {
-    //     result.ascanBlockPreprocessed = ascanBlock.ascanBlock;
-    // }
-//printTime();
-//std::cout<<"end"<<std::endl;    
-    return result;
-}                                                         
+}                                                        
--- a/src/common/dataBlockCreation/getAScanBlockPreprocessed.h
+++ b/src/common/dataBlockCreation/getAScanBlockPreprocessed.h
@@ -26,23 +26,6 @@ namespace Recon
    AscanBlockPreprocessed getAscanBlockPreprocessed(Parser* aParser, const Aurora::Matrix& aMp, const Aurora::Matrix& aSl, const Aurora::Matrix& aSn, 
                                                     const Aurora::Matrix& aRl, const Aurora::Matrix& aRn, GeometryInfo& aGeom, const MeasurementInfo& aMeasInfo,
                                                     bool aApplyFilter, bool aTransReco);
-
-    struct AscanBlockPreprocessedCuda
-    {
-        Aurora::CudaMatrix ascanBlockPreprocessed;
-        Aurora::Matrix mpBlock;
-        Aurora::Matrix slBlock;
-        Aurora::Matrix snBlock;
-        Aurora::Matrix rlBlock;
-        Aurora::Matrix rnBlock;
-        Aurora::Matrix senderPositionBlock;
-        Aurora::Matrix receiverPositionBlock;
-        Aurora::CudaMatrix gainBlock;
-    };
-
-    AscanBlockPreprocessedCuda getAscanBlockPreprocessedCuda(Parser* aParser, const Aurora::Matrix& aMp, const Aurora::Matrix& aSl, const Aurora::Matrix& aSn, 
-                                                         const Aurora::Matrix& aRl, const Aurora::Matrix& aRn, GeometryInfo& aGeom, const MeasurementInfo& aMeasInfo,
-                                                         bool aApplyFilter, bool aTransReco);
 }


--- a/src/common/fileHelper.h
+++ b/src/common/fileHelper.h
@@ -6,7 +6,7 @@
 namespace Recon
 {
    const std::string DEFAULT_CONFIG_PATH = "/home/UR/ConfigFiles/";
-    const std::string DEFAULT_OUTPUT_PATH = "/home/UR/ReconResult/";
+    const std::string DEFAULT_OUTPUT_PATH = "./";
    const std::string DEFAULT_OUTPUT_FILENAME = "sun.mat";

    std::string getPath(const std::string &aFullPath);
--- a/src/common/getMeasurementMetaData.cpp
+++ b/src/common/getMeasurementMetaData.cpp
@@ -14,6 +14,7 @@
 #include "Function1D.h"
 #include "Function2D.h"
 #include "AuroraDefs.h"
+#include "log/log.h"
 #include "mkl.h"

 #include "../config/config.h"
@@ -275,7 +276,7 @@ CEInfo Recon::getCEInfo(Parser* aParser, const MeasurementInfo aInfo)
        short* ceMeasuredData = ceMeasuredPointer.get();
        if(aParser->getMetaData().getDataType() == "float16")
        {
-            ce = convertfp16tofloat(ceMeasuredPointer.get(), ceMeasuredLength, ceMeasuredPointer.getSize());
+            ce = convertfp16tofloat(ceMeasuredPointer.get(), ceMeasuredLength, ceMeasuredPointer.getSize()); 
        }
        else
        {
--- a/src/common/preprocessAscanBlock.cpp
+++ b/src/common/preprocessAscanBlock.cpp
@@ -21,14 +21,4 @@ Aurora::Matrix Recon::preprocessAscanBlock(const Aurora::Matrix& aAscans, const
    // end

    return result;
-}
-
-Aurora::CudaMatrix Recon::preprocessAscanBlockCuda(const Aurora::CudaMatrix& aAscans, const MeasurementInfo& aMeasInfo)
-{
-    if(aMeasInfo.ascanDataType == "float16")
-    {
-        return Aurora::convertfp16tofloatCuda(aAscans, aAscans.getDimSize(0), aAscans.getDimSize(1));
-    }
-
-    return aAscans;
 }
--- a/src/common/preprocessAscanBlock.h
+++ b/src/common/preprocessAscanBlock.h
@@ -8,7 +8,6 @@ namespace Recon
 {
    Aurora::Matrix preprocessAscanBlock(const Aurora::Matrix& aAscans, const MeasurementInfo& aMeasInfo);

-    Aurora::CudaMatrix preprocessAscanBlockCuda(const Aurora::CudaMatrix& aAscans, const MeasurementInfo& aMeasInfo);
 }

 #endif
--- a/src/config/config.cpp
+++ b/src/config/config.cpp
@@ -307,11 +307,11 @@ namespace Recon
                    nlohmann::json transmissionCorrection = params.at("transmissionCorrection");
                    if(transmissionCorrection.contains("soundSpeedCorrection"))
                    {
-                        reflectParams::soundSpeedCorrection = transmissionCorrection.at("soundSpeedCorrection").get<int>();
+                        reflectParams::soundSpeedCorrection = 0;
                    }
                    if(transmissionCorrection.contains("attenuationCorrection"))
                    {
-                        reflectParams::attenuationCorrection = transmissionCorrection.at("attenuationCorrection").get<int>();
+                        reflectParams::attenuationCorrection = 0;
                    }
                    if(transmissionCorrection.contains("saveTransmInReflCoords"))
                    {
@@ -679,7 +679,7 @@ namespace Recon
        reflectParams::windowLength = 10;
        reflectParams::numThreads = 30;
        reflectParams::expectedUSSpeedRange = Aurora::Matrix::fromRawData(new float[2] {1420,1600}, 1, 2);
-        //reflectParams.transmissionCorrection
+        // reflectParams.transmissionCorrection
        reflectParams::soundSpeedCorrection = 0;
        reflectParams::attenuationCorrection = 0;
        reflectParams::saveTransmInReflCoords = 1;
--- a/src/main.cxx
+++ b/src/main.cxx
@@ -12,7 +12,6 @@
 #include "startReconstructions.h"

 #include "transmissionReconstruction/detection/detection.h"
-#include "transmissionReconstruction/detection/detection.cuh"
 #include "startReconstructions.h"
 #include "log/log.h"
 /* 0 is data path.
--- a/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp
+++ b/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp
@@ -7,6 +7,7 @@
 #include "Matrix.h"
 #include "common/dataBlockCreation/removeDataFromArrays.h"
 #include "config/config.h"
+#include "log/log.h"
 #include <algorithm>
 #include <cmath>
 #include <cstddef>
@@ -14,6 +15,7 @@
 #include <iostream>
 #include <memory>
 #include <mkl_cblas.h>
+#include <mkl_service.h>
 #include <mkl_vml_functions.h>
 #include <vector>
 namespace Recon {
@@ -120,6 +122,7 @@ namespace Recon {
    const Aurora::Matrix &blockedReceiverPosition,
    const Aurora::Matrix &blockedGain, const Aurora::Matrix &blockedChannels,
    MeasurementInfo &info, const PreComputes &preComputes) {
+        RECON_INFO("============================>performSignalProcessing start!");
        std::vector<Aurora::Matrix> result;
        auto t = size(blockedAScans);
        size_t nSamples = (int)t[0];
@@ -132,9 +135,13 @@ namespace Recon {
        }
        int winLength = 100; 
        auto ascanMapValue = Aurora::zeros(1,blockedAScans.getDimSize(1),1);
-        #pragma omp parallel for num_threads(32)
+        RECON_INFO("performSignalProcessing step 1 done!");
+
+        // #pragma omp parallel for num_threads(32)
+        #pragma omp parallel for 
        for (size_t i = 0; i < blockedAScans.getDimSize(1); i++)
        {
+            // mkl_set_num_threads_local(1);
            auto _blockedAScans = getPartMatrixColRef(blockedAScans,i);
            _blockedAScans = _blockedAScans/blockedGain[i];
            if (reflectParams::removeDCOffset == 1)
@@ -176,11 +183,15 @@ namespace Recon {
            }
            blockedAScans(Aurora::$,i)= _blockedAScans;
        }
+        RECON_INFO("performSignalProcessing loop 1 done!");
+
        auto fx_all = fft(blockedAScans);
+        RECON_INFO("performSignalProcessing fft 1 done!");

        if (reflectParams::useCorrelation && reflectParams::matchedFilterCeAScan)
        {
-            #pragma omp parallel for num_threads(32)
+            // #pragma omp parallel for num_threads(32)
+            #pragma omp parallel for 
            for (size_t i = 0; i < blockedAScans.getDimSize(1); i++)
            {
                Aurora::Matrix fx = getPartMatrixColRef(fx_all,i);
@@ -214,26 +225,41 @@ namespace Recon {
                Aurora::free(value2);
                fx_all(Aurora::$,i) = complex;
            }
+            RECON_INFO("performSignalProcessing loop 2 done!");
+
            blockedAScans = Aurora::real(ifft(fx_all));
+            RECON_INFO("performSignalProcessing ifft 2 done!");

            blockedAScans.setBlockValue(0, 0, winLength-1, 0);
            blockedAScans.setBlockValue(0, blockedAScans.getDimSize(0)-winLength, blockedAScans.getDimSize(0)-1, 0);
            blockedAScans = fft(blockedAScans);
+            RECON_INFO("performSignalProcessing fft 2 done!");
+
        }
        else{
            blockedAScans = fx_all;
        }
+        RECON_INFO("->offsetSignalFourier start!");
+
        auto exponent = offsetSignalFourier(preComputes, nSamples);
        // exponent.forceReshape(1, exponent.getDataSize(), 1);
-        #pragma omp parallel for num_threads(32)
+        RECON_INFO("->offsetSignalFourier finish!");
+
+        // #pragma omp parallel for num_threads(32)
+        #pragma omp parallel for 
        for (size_t i = 0; i < blockedAScans.getDimSize(1); i++)
        {
            Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i);
            blockedAScans(Aurora::$,i) =_blockedAScans*exponent;
        }
+        RECON_INFO("performSignalProcessing loop 3 done!");
+
        blockedAScans = real(ifft(blockedAScans));
+        RECON_INFO("performSignalProcessing ifft 3 done!");
+
        if (reflectParams::removeTransmissionSignal == 1) {
-            #pragma omp parallel for num_threads(32)
+            // #pragma omp parallel for num_threads(32)
+            #pragma omp parallel for 
            for (size_t i = 0; i < blockedAScans.getDimSize(1); i++)
            {
                Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i);
@@ -333,13 +359,18 @@ namespace Recon {
                    // blockedAScans(Aurora::$,i) =_blockedAScans;
                }
            }
+            RECON_INFO("performSignalProcessing loop 4 done!");
+
        }
        // auto __blockedAScans = fft(blockedAScans(Aurora::$,0).toMatrix());
        blockedAScans = fft(blockedAScans);
+        RECON_INFO("performSignalProcessing fft 4 done!");
+
        if(reflectParams::useOptPulse==1){
            auto n = nSamples;
            auto nHalf = round((float)n/2);
-            #pragma omp parallel for num_threads(32)
+            // #pragma omp parallel for num_threads(32)
+            #pragma omp parallel for 
            for (size_t i = 0; i < blockedAScans.getDimSize(1); i++)
            {
                Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i);
@@ -348,13 +379,18 @@ namespace Recon {
                _blockedAScans.setBlockComplexValue(0,nHalf,_blockedAScans.getDimSize(0)-1, std::complex<float>{0,0});
                _blockedAScans.setBlock(0,1,nHalf-1,_blockedAScans.block(0,1,nHalf-1)*2);
            }
+            RECON_INFO("performSignalProcessing loop 5 done!");
        
            blockedAScans = ifft(blockedAScans, n);
+            RECON_INFO("performSignalProcessing ifft 5 done!");
+
        
            blockedAScans = abs(blockedAScans);

            auto help_all = blockedAScans.deepCopy();
-            #pragma omp parallel for num_threads(32)
+
+            // #pragma omp parallel for num_threads(32)
+            #pragma omp parallel for 
            for (size_t i = 0; i < blockedAScans.getDimSize(1); i++)
            {
                Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i);
@@ -395,28 +431,40 @@ namespace Recon {
                }
                help_all(Aurora::$,i) = help;
            }
+            RECON_INFO("performSignalProcessing loop 6 done!");

            help_all = fft(help_all);
+            RECON_INFO("performSignalProcessing fft 6 done!");

-            #pragma omp parallel for num_threads(32)
+
+            // #pragma omp parallel for num_threads(32)
+            #pragma omp parallel for 
            for (size_t i = 0; i < help_all.getDimSize(1); i++)
            {
+                // mkl_set_num_threads_local(1);
                Aurora::Matrix help = getPartMatrixColRef(help_all,i);
                help_all(Aurora::$,i) = help*(preComputes.sincPeak_ft);
            }
+            RECON_INFO("performSignalProcessing loop 7 done!");

            help_all = real(ifft(help_all));
+            RECON_INFO("performSignalProcessing ifft 7 done!");

-            #pragma omp parallel for num_threads(32)
+            // #pragma omp parallel for num_threads(32)
+            #pragma omp parallel for 
            for (size_t i = 0; i < help_all.getDimSize(1); i++)
            {
                Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i);
                Aurora::Matrix help = getPartMatrixColRef(help_all,i);
                Aurora::compareSet(_blockedAScans,0,help,Aurora::NL);
            }
+            RECON_INFO("performSignalProcessing loop 8 done!");
+
        }
        else{
             blockedAScans = real(ifft(blockedAScans));
+            RECON_INFO("performSignalProcessing ifft 8.1 done!");
+
        }
        result.push_back(blockedAScans);
        
@@ -427,6 +475,7 @@ namespace Recon {
        }
        
        result.push_back(ascanMapValue);
+        RECON_INFO("============================>performSignalProcessing finish!");
        return result;
    }

--- a/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp
+++ b/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp
@@ -4,9 +4,6 @@
 #include "Function2D.h"
 #include "Function3D.h"
 #include "Matrix.h"
-#include "SAFTStructs.h"
-// #include "SAFT_ATT.h"
-#include "SAFT_TOFI.h"
 #include "common/dataBlockCreation/removeDataFromArrays.h"
 #include "config/config.h"

@@ -19,314 +16,4 @@
 #include <vector>

 namespace Recon {
-
-    int compareColumn(const Aurora::Matrix &aMatrix, int colCount,
-                    const Aurora::Matrix &aOther, int col2) {
-
-    for (size_t j = 0; j < colCount; j++) {
-        size_t i = 0;
-        for (; i < aMatrix.getDimSize(0); i++) {
-        if (aMatrix[j * aMatrix.getDimSize(0) + 1] !=
-            aOther[col2 * aOther.getDimSize(0) + i])
-            break;
-        }
-        if (i == aMatrix.getDimSize(0))
-        return j;
-    }
-    return -1;
-    }
-
-
-    Aurora::Matrix callSAFT(const Aurora::Matrix &AScans,
-                            const Aurora::Matrix &receiverIndex,
-                            const Aurora::Matrix &senderIndex,
-                            const Aurora::Matrix &receiverPosGeom,
-                            const Aurora::Matrix &senderPosGeom,
-                            int SAFT_mode, float TimeInterval,
-                            const TransRecos& transRecos,
-                            const Aurora::Matrix &Env) 
-    {
-        if (reflectParams::useAscanIndex == 0) {
-            if (reflectParams::soundSpeedCorrection == 0) {
-                std::cerr
-                    << "reconstruction without sound speed correction not allowed!!!!"
-                    << std::endl;
-                return Aurora::Matrix();
-            }
-        }
-        std::vector<Matrix_t> params;
-        Matrix_t AScan{nullptr,
-                        (size_t)AScans.getDims(),
-                        {(size_t)AScans.getDimSize(0), (size_t)AScans.getDimSize(1),
-                            (size_t)AScans.getDimSize(2)},
-                        AScans.getDataSize()};
-        float *ASdata = new float[AScans.getDataSize()]{0};
-        std::copy(AScans.getData(), AScans.getData() + AScans.getDataSize(),
-                ASdata);
-        AScan.Data = ASdata;
-        params.push_back(AScan);
-
-        Matrix_t imageStartpoint{nullptr, 2, {1, 3}, 3};
-        imageStartpoint.Data =
-            new float[3]{(float)reflectParams::imageStartpoint[0],
-                        (float)reflectParams::imageStartpoint[1],
-                        (float)reflectParams::imageStartpoint[2]};
-        params.push_back(imageStartpoint);
-
-        Matrix_t receiverIdx{nullptr,
-                            2,
-                            { 1,receiverIndex.getDataSize(), 1},
-                            receiverIndex.getDataSize()};
-        unsigned short *temp1 = new unsigned short[receiverIndex.getDataSize()];
-        std::copy(receiverIndex.getData(),
-                receiverIndex.getData() + receiverIndex.getDataSize(), temp1);
-        receiverIdx.Data = temp1;
-        params.push_back(receiverIdx);
-
-        Matrix_t senderIdx{nullptr,
-                        2,
-                        {1,senderIndex.getDataSize(), 1},
-                        senderIndex.getDataSize()};
-        auto temp2 = new unsigned short[senderIndex.getDataSize()];
-        std::copy(senderIndex.getData(),
-                senderIndex.getData() + senderIndex.getDataSize(), temp2);
-        senderIdx.Data = temp2;
-        params.push_back(senderIdx);
-
-        Matrix_t receiverGeom{nullptr,
-                            (size_t)receiverPosGeom.getDims(),
-                            {(size_t)receiverPosGeom.getDimSize(0),
-                            (size_t)receiverPosGeom.getDimSize(1),
-                            (size_t)receiverPosGeom.getDimSize(2)},
-                            receiverPosGeom.getDataSize()};
-        float *fdata = new float[receiverPosGeom.getDataSize()]{0};
-        std::copy(receiverPosGeom.getData(),
-                receiverPosGeom.getData() + receiverPosGeom.getDataSize(), fdata);
-        receiverGeom.Data = fdata;
-        params.push_back(receiverGeom);
-
-        Matrix_t senderGeom{nullptr,
-                            (size_t)senderPosGeom.getDims(),
-                            {(size_t)senderPosGeom.getDimSize(0),
-                            (size_t)senderPosGeom.getDimSize(1),
-                            (size_t)senderPosGeom.getDimSize(2)},
-                            senderPosGeom.getDataSize()};
-        auto fdata1 = new float[senderPosGeom.getDataSize()]{0};
-        std::copy(senderPosGeom.getData(),
-                senderPosGeom.getData() + senderPosGeom.getDataSize(), fdata1);
-        senderGeom.Data = fdata1;
-        params.push_back(senderGeom);
-
-        Matrix_t _SAFT_mode{&SAFT_mode, 1, {1, 1, 1}, 1};
-        params.push_back(_SAFT_mode);
-
-        int saftVariant_v[6]{(int)reflectParams::saftVariant[0],(int)reflectParams::saftVariant[1],(int)reflectParams::saftVariant[2],
-                            (int)reflectParams::saftVariant[3],(int)reflectParams::saftVariant[4],(int)reflectParams::saftVariant[5]};
-        Matrix_t saftVariant{saftVariant_v,2,{1,6,1},6};
-        params.push_back(saftVariant);
-
-        Matrix_t SpeedMap3D{nullptr,
-                            (size_t)transRecos.speedMap3D.getDims(),
-                            {(size_t)transRecos.speedMap3D.getDimSize(0),
-                            (size_t)transRecos.speedMap3D.getDimSize(1),
-                            (size_t)transRecos.speedMap3D.getDimSize(2)},
-                            transRecos.speedMap3D.getDataSize()};
-        auto fdata2 = new float[transRecos.speedMap3D.getDataSize()]{0};
-        std::copy(transRecos.speedMap3D.getData(),
-                transRecos.speedMap3D.getData() + transRecos.speedMap3D.getDataSize(), fdata2);
-        SpeedMap3D.Data = fdata2;
-        params.push_back(SpeedMap3D);
-
-        Matrix_t BeginTransMap{nullptr,
-                            (size_t)transRecos.beginTransMap.getDims(),
-                            {(size_t)transRecos.beginTransMap.getDimSize(0),
-                            (size_t)transRecos.beginTransMap.getDimSize(1),
-                            (size_t)transRecos.beginTransMap.getDimSize(2)},
-                            transRecos.beginTransMap.getDataSize()};
-        auto fdata3 = new float[transRecos.beginTransMap.getDataSize()]{0};
-        std::copy(transRecos.beginTransMap.getData(),
-                transRecos.beginTransMap.getData() + transRecos.beginTransMap.getDataSize(), fdata3);
-        BeginTransMap.Data = fdata3;
-        params.push_back(BeginTransMap);
-
-        Matrix_t DeltaTransMap{nullptr,
-                            1,
-                            1,
-                            1,
-                            1,
-                            1};
-        auto fdata4 = new float[1]{0};
-        fdata4[0] = transRecos.deltaTransMap;
-        DeltaTransMap.Data = fdata4;
-        params.push_back(DeltaTransMap);
-
-        Matrix_t AttenuationMap3D{nullptr,
-                            (size_t)transRecos.attenuationMap3D.getDims(),
-                            {(size_t)transRecos.attenuationMap3D.getDimSize(0),
-                            (size_t)transRecos.attenuationMap3D.getDimSize(1),
-                            (size_t)transRecos.attenuationMap3D.getDimSize(2)},
-                            transRecos.attenuationMap3D.getDataSize()};
-        auto fdata5 = new float[transRecos.attenuationMap3D.getDataSize()]{0};
-        std::copy(transRecos.attenuationMap3D.getData(),
-                transRecos.attenuationMap3D.getData() + transRecos.attenuationMap3D.getDataSize(), fdata5);
-        AttenuationMap3D.Data = fdata5;
-        params.push_back(AttenuationMap3D);
-
-        float imageResolution_v = (float)reflectParams::imageResolution;
-        Matrix_t imageResolution{&imageResolution_v, 1, {1, 1, 1}, 1};
-        params.push_back(imageResolution);
-        
-        float TimeInterval_v = (float)TimeInterval;
-        Matrix_t _TimeInterval{&TimeInterval_v, 1, {1, 1, 1}, 1};
-        params.push_back(_TimeInterval);
-
-        int imageXYZ_v[3]{(int)reflectParams::imageXYZ[0],(int)reflectParams::imageXYZ[1],(int)reflectParams::imageXYZ[2]};
-        Matrix_t imageXYZ{imageXYZ_v,2,{1,3,1},3};
-        params.push_back(imageXYZ);
-
-        Matrix_t _Env{nullptr,
-                    (size_t)Env.getDims(),
-                    {(size_t)Env.getDimSize(0), (size_t)Env.getDimSize(1),
-                        (size_t)Env.getDimSize(2)},
-                    Env.getDataSize()};
-        _Env.Data = Env.getData();
-        params.push_back(_Env);
-
-        int blockXYZ_v[3]{(int)reflectParams::blockDimXYZ[0],(int)reflectParams::blockDimXYZ[1],(int)reflectParams::blockDimXYZ[2]};
-        Matrix_t blockXYZ{blockXYZ_v,2,{1,3,1},3};
-        params.push_back(blockXYZ);
-
-        Matrix_t gpus{nullptr,
-                    2,
-                    {1, (size_t)reconParams::gpuSelectionList.getDataSize(),
-                        1},
-                    reconParams::gpuSelectionList.getDataSize()};
-        int *gdata = new int[reconParams::gpuSelectionList.getDataSize()]{0};
-        std::copy(reconParams::gpuSelectionList.getData(), reconParams::gpuSelectionList.getData() + reconParams::gpuSelectionList.getDataSize(),
-                gdata);
-        gpus.Data = gdata;
-        params.push_back(gpus);
-
-        if (reflectParams::useAscanIndex == 0){
-            Matrix_t medianWindowSize{&reflectParams::medianWindowSize,1,{1,1,1},1};
-            params.push_back(medianWindowSize);
-        }
-
-        float other[2]={(float)reflectParams::debugMode, (float)reflectParams::attenuationCorrectionLimit} ;
-        Matrix_t otherP{other,2,{1,2,1},2};
-        params.push_back(otherP);
-        if (reflectParams::useAscanIndex == 0) {
-            // auto result =  SAFT_ATT(params);
-            // delete [] ASdata;
-            // delete [] (float*)imageStartpoint.Data;
-            // delete [] temp1;
-            // delete [] temp2;
-            // delete [] fdata;
-            // delete [] fdata1;
-            // delete [] fdata2;
-            // delete [] fdata3;
-            // delete [] fdata4;
-            // delete [] fdata5;
-            // delete [] gdata;
-            // return Aurora::Matrix::fromRawData((float *)result.Data, result.Dims[0],
-            //                                 result.Dims[1], result.Dims[2]);
-
-            return Aurora::Matrix();
-        } else {
-            auto result = SAFT_TOFI(params);
-            delete [] ASdata;
-            delete [] (float*)imageStartpoint.Data;
-            delete [] temp1;
-            delete [] temp2;
-            delete [] fdata;
-            delete [] fdata1;
-            delete [] fdata2;
-            delete [] fdata3;
-            delete [] fdata4;
-            delete [] fdata5;
-            delete [] gdata;
-            return Aurora::Matrix::fromRawData((float *)result.Data, result.Dims[0],
-                                            result.Dims[1], result.Dims[2]);
-        }
-    }
-
-    Aurora::Matrix
-    recontructSAFT(const Aurora::Matrix &AScans, const Aurora::Matrix &senderPos,
-                const Aurora::Matrix &receiverPos, const Aurora::Matrix &mpIndex,
-                int SAFT_mode, TransRecos &transRecos, Aurora::Matrix &Env) {
-        float TimeInterval = 1.0 / (reflectParams::aScanReconstructionFrequency);
-        std::vector<int> motorPosAvailable;
-        std::unique_copy(mpIndex.getData(), mpIndex.getData() + mpIndex.getDataSize(),
-                         std::back_inserter(motorPosAvailable));
-        for (auto mp : motorPosAvailable) {
-            auto mpIdxs = mpIndex == mp;
-            size_t countMp = Aurora::sum(mpIdxs).getScalar();
-
-            auto senderPos_s = Aurora::transpose(removeDataFromArrays(senderPos, mpIdxs));
-            auto receiverPos_s = Aurora::transpose(removeDataFromArrays(receiverPos, mpIdxs));
-            Aurora::Matrix idsSender ;
-            auto senderPosGeom = Aurora::transpose(Aurora::uniqueByRows(senderPos_s, idsSender));
-            idsSender = Aurora::transpose(idsSender);
-            Aurora::Matrix idsReceiver ;
-            auto receiverPosGeom = Aurora::transpose(Aurora::uniqueByRows(receiverPos_s, idsReceiver));
-            idsReceiver = Aurora::transpose(idsReceiver);
-
-            auto numUsedData = countMp;
-            int count = numUsedData / reflectParams::blockSizeReco;
-            auto blockIdxs = Aurora::zeros(1, count + 2);
-            for (size_t i = 0; i < count + 1; i++) {
-                blockIdxs[i] = i * reflectParams::blockSizeReco;
-            }
-            blockIdxs[count + 1] = numUsedData;
-            for (size_t i = 0; i < count + 1; i++) {
-                size_t length = blockIdxs[i + 1] - blockIdxs[i];
-                size_t begin = blockIdxs[i];
-                size_t end = blockIdxs[i + 1];
-                
-                auto receiverIndex = Aurora::zeros(1, length);
-                auto senderIndex = Aurora::zeros(1, length);
-                auto mpIdxs_d = Aurora::zeros(1,mpIdxs.getDataSize());
-                for (size_t j = begin, k = 0; j < end; j++, k++) {
-                    mpIdxs_d[j] =mpIdxs[j];
-                    receiverIndex[k] = idsReceiver[j]+1;
-                    senderIndex[k] = idsSender[j]+1;
-                }
-                auto _AScans = removeDataFromArrays(AScans, mpIdxs_d);
-                Env = callSAFT(_AScans, receiverIndex, senderIndex, receiverPosGeom,
-                            senderPosGeom, SAFT_mode, TimeInterval, transRecos, Env);
-            }
-        }
-        return Env;
-    }
-
-    //TODO:untested
-    Aurora::Matrix polyvaln(PolyModel polymodel, const Aurora::Matrix &indepvar)
-    {
-        auto np = Aurora::size(indepvar);
-        int n = (int)np[0];
-        int p = (int)np[1];
-        Aurora::Matrix _indepvar = indepvar;
-        if (n == 1 && polymodel.ModelTerms.getDimSize(1) == 1){
-            _indepvar = Aurora::transpose(indepvar);
-            np = Aurora::size(indepvar);
-            n = (int)np[0];
-            p = (int)np[1];
-        }
-        else if(polymodel.ModelTerms.getDimSize(1) != p){
-            std::cerr<<"Size of indepvar array and this model are inconsistent."<<std::endl;
-            return Aurora::Matrix();
-        }
-        size_t nt = Aurora::size(polymodel.ModelTerms,1);
-        auto ypred = Aurora::zeros(n,1);
-        for (size_t i = 0; i < nt; i++)
-        {
-            auto t = Aurora::ones(n,1);
-            for (size_t j = 0; j < p; j++)
-            {
-                t = (t*(_indepvar(Aurora::$,j).toMatrix()^polymodel.Coefficients(i,j).toMatrix().getScalar()));
-            }
-            ypred = ypred+t*polymodel.Coefficients[i];
-        }
-        return ypred;
-    }
 } // namespace Recon
--- a/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.h
+++ b/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.h
@@ -3,18 +3,6 @@
 #include "Matrix.h"
 #include "src/reflectionReconstruction/preprocessData/preprocessTransmissionReconstructionForReflection.h"
 namespace Recon {
-
-struct PolyModel{
-    Aurora::Matrix ModelTerms;
-    Aurora::Matrix Coefficients;
-};
-
-Aurora::Matrix
-recontructSAFT(const Aurora::Matrix &AScans, const Aurora::Matrix &senderPos,
-               const Aurora::Matrix &receiverPos, const Aurora::Matrix &mpIndex,
-               int SAFT_mode, TransRecos &transRecos, Aurora::Matrix &Env);
-
-Aurora::Matrix polyvaln(PolyModel polymodel, const Aurora::Matrix &indepvar);
 }

 #endif // __RECONSTRUCTIONSAFT_H__
--- a/src/reflectionReconstruction/startReflectionReconstruction.cpp
+++ b/src/reflectionReconstruction/startReflectionReconstruction.cpp
@@ -16,10 +16,9 @@
 #include "config/config.h"
 #include "log/log.h"

-#include "CudaEnvInit.h"
-
 #include <cstdio>
 #include <iostream>
+#include <ostream>
 #include <queue>
 #include <thread>

@@ -39,20 +38,17 @@ void producerThread( Parser* aParser, const Aurora::Matrix& aMotorPos,
                     const Aurora::Matrix& aRlList, const Aurora::Matrix& aRnList,
                     GeometryInfo aGeom, MeasurementInfo aExpInfo, PreComputes aPreComputes)
 {
-printf("Reflection reconstruction is carried out.");
+    // printf("Reflection reconstruction is carried out.");

-    printf("Preperations for reconstructions.");
+    // printf("Preperations for reconstructions.");

-    printf("   - reset GPUs");
+    // printf("   - reset GPUs");
    for (size_t i = 0; i < reflectParams::gpuSelectionList.getDataSize(); i++)
    {
        std::string msg;
-        if (!resetGPUDevice((int)reflectParams::gpuSelectionList[i],msg))
-        {
-            std::cerr<<msg<<std::endl;
-        }
+
    }
-    printf("   - optimal pulse");
+    // printf("   - optimal pulse");
    if(reflectParams::useOptPulse==1 && reflectParams::runReflectionReco)
    {
        aPreComputes.sincPeak_ft = determineOptimalPulse(aPreComputes.timeInterval, aExpInfo.expectedAScanLength);
@@ -62,17 +58,25 @@ printf("Reflection reconstruction is carried out.");
                      aSnList.getDataSize() * aRlList.getDataSize() *
                      aRnList.getDataSize();
    int numTakenScans = 0,numProcessedScans = 0,numPossibleScans = 0;
+    int index=1;
    for(int i=0; i<aMotorPos.getDataSize(); ++i)
    {
-        for(int j=0; j<aSlList.getDataSize() / transParams::senderTASSize; ++j)
+        for(int j=0; j</*aSlList.getDataSize() / transParams::senderTASSize*/1; ++j)
        {
-            for(int k=0; k<aSnList.getDataSize() / transParams::senderElementSize; ++k)
+            for(int k=0; k</*aSnList.getDataSize() / transParams::senderElementSize*/1; ++k)
            {
+    // for(int i=0; i<1; ++i)
+    // {
+    //     for(int j=0; j<1; ++j)
+    //     {
+    //         for(int k=0; k<3; ++k)
+    //         {
                Matrix mp = aMotorPos(i).toMatrix();
                Matrix sl = aSlList.block(0, transParams::senderTASSize*j, transParams::senderTASSize*j+transParams::senderTASSize - 1);
                Matrix sn = aSnList.block(0, transParams::senderElementSize*k, transParams::senderElementSize*k+transParams::senderElementSize - 1);
+                RECON_INFO("getAscanBlockPreprocessed start");
                auto blockData = getAscanBlockPreprocessed(aParser, mp, sl, sn, aRlList, aRnList, aGeom, aExpInfo, true, false);
-
+                RECON_INFO("getAscanBlockPreprocessed end");
                float* channelListSizeData = Aurora::malloc(2);
                channelListSizeData[0] = channelList.getDimSize(0);
                channelListSizeData[1] = channelList.getDimSize(1);
@@ -85,15 +89,18 @@ printf("Reflection reconstruction is carried out.");
                    channelBlockData[i] = channelList[ind[i] - 1];
                }
                Matrix channelBlock = Matrix::New(channelBlockData, 1, channelBlockSize);
-                                auto preprocessData = preprocessAScanBlockForReflection(blockData.ascanBlockPreprocessed, blockData.mpBlock, blockData.slBlock,
+                RECON_INFO("preprocessAScanBlockForReflection start");
+                auto preprocessData = preprocessAScanBlockForReflection(blockData.ascanBlockPreprocessed, blockData.mpBlock, blockData.slBlock,
                                      blockData.snBlock, blockData.rlBlock, blockData.rnBlock, blockData.senderPositionBlock,
                                      blockData.receiverPositionBlock, blockData.gainBlock, channelBlock, aExpInfo, aPreComputes);
-
-                std::unique_lock<std::mutex> lock(PRODUCER_MUTEX); 
-                PRODUCER_BLOCKDATAS.push(blockData);
-                PRODUCER_PROCESSDATAS.push(preprocessData);
-                lock.unlock();
-                PRODUCER_CONDITION.notify_one();
+                RECON_INFO("preprocessAScanBlockForReflection end");
+                RECON_INFO("Reflect reconstructing...... {0}",index);
+                index++;
+                // std::unique_lock<std::mutex> lock(PRODUCER_MUTEX); 
+                // PRODUCER_BLOCKDATAS.push(blockData);
+                // PRODUCER_PROCESSDATAS.push(preprocessData);
+                // lock.unlock();
+                // PRODUCER_CONDITION.notify_one();
            }
        }
    }
@@ -108,41 +115,41 @@ Aurora::Matrix Recon::startReflectionReconstruction( Parser* aParser, int aSAFT_
    for (size_t i = 0; i < reflectParams::gpuSelectionList.getDataSize(); i++)
    {
        std::string msg;
-        if (!resetGPUDevice((int)reflectParams::gpuSelectionList[i],msg))
-        {
-            std::cerr<<msg<<std::endl;
-        }
    }
-    std::thread thread = std::thread(producerThread, aParser, aMotorPos, aSlList, aSnList, aRlList, aRnList, aGeom, aExpInfo, aPreComputes);
+    // std::thread thread = std::thread(producerThread, aParser, aMotorPos, aSlList, aSnList, aRlList, aRnList, aGeom, aExpInfo, aPreComputes);
+    int index=1;
+    // for(int i=0; i<aMotorPos.getDataSize(); ++i)
+    // {
+    //     for(int j=0; j<aSlList.getDataSize() / transParams::senderTASSize; ++j)
+    //     {
+    //         for(int k=0; k<aSnList.getDataSize() / transParams::senderElementSize; ++k)
+    //         {
+                producerThread(aParser, aMotorPos, aSlList, aSnList, aRlList, aRnList, aGeom, aExpInfo, aPreComputes);
+    // for(int i=0; i<1; ++i)
+    // {
+    //     for(int j=0; j<1; ++j)
+    //     {
+    //         for(int k=0; k<3; ++k)
+    //         {
+                // std::unique_lock<std::mutex> lock(PRODUCER_MUTEX); 
+                // PRODUCER_CONDITION.wait(lock, []{return !PRODUCER_PROCESSDATAS.empty() && !PRODUCER_BLOCKDATAS.empty();});
+                // lock.unlock();
+                // Env = recontructSAFT(removeDataFromArrays(PRODUCER_PROCESSDATAS.front().AscanBlock, PRODUCER_PROCESSDATAS.front().usedData),
+                //                      removeDataFromArrays(PRODUCER_BLOCKDATAS.front().senderPositionBlock, PRODUCER_PROCESSDATAS.front().usedData),
+                //                      removeDataFromArrays(PRODUCER_BLOCKDATAS.front().receiverPositionBlock, PRODUCER_PROCESSDATAS.front().usedData),
+                //                      removeDataFromArrays(PRODUCER_BLOCKDATAS.front().mpBlock, PRODUCER_PROCESSDATAS.front().usedData),
+                //                      aSAFT_mode, aTransRecos, Env);
+                // lock.lock();
+                // PRODUCER_PROCESSDATAS.pop();
+                // PRODUCER_BLOCKDATAS.pop();
+                // lock.unlock();

-    Matrix Env =  Aurora::zeros((int)reflectParams::imageXYZ[0],(int)reflectParams::imageXYZ[1],(int)reflectParams::imageXYZ[2]);
+    //         }
+    //         Recon::notifyProgress(25+73*((j*i)/(aMotorPos.getDataSize() * aSlList.getDataSize())));
+    //     }
+    // }

-    for(int i=0; i<aMotorPos.getDataSize(); ++i)
-    {
-        //#pragma omp parallel for num_threads(24)
-        for(int j=0; j<aSlList.getDataSize() / transParams::senderTASSize; ++j)
-        {
-            for(int k=0; k<aSnList.getDataSize() / transParams::senderElementSize; ++k)
-            {
-                std::unique_lock<std::mutex> lock(PRODUCER_MUTEX); 
-                PRODUCER_CONDITION.wait(lock, []{return !PRODUCER_PROCESSDATAS.empty() && !PRODUCER_BLOCKDATAS.empty();});
-                lock.unlock();
-                Env = recontructSAFT(removeDataFromArrays(PRODUCER_PROCESSDATAS.front().AscanBlock, PRODUCER_PROCESSDATAS.front().usedData),
-                                     removeDataFromArrays(PRODUCER_BLOCKDATAS.front().senderPositionBlock, PRODUCER_PROCESSDATAS.front().usedData),
-                                     removeDataFromArrays(PRODUCER_BLOCKDATAS.front().receiverPositionBlock, PRODUCER_PROCESSDATAS.front().usedData),
-                                     removeDataFromArrays(PRODUCER_BLOCKDATAS.front().mpBlock, PRODUCER_PROCESSDATAS.front().usedData),
-                                     aSAFT_mode, aTransRecos, Env);
-                lock.lock();
-                PRODUCER_PROCESSDATAS.pop();
-                PRODUCER_BLOCKDATAS.pop();
-                lock.unlock();
-                std::cout<<Env[0]<<"-" << Env[1] <<"-" << Env[2] <<"-" << Env[3]<<std::endl;
-                RECON_INFO("Reflection Reconstructon: " + std::to_string(j));
-            }
-            Recon::notifyProgress(25+73*((j*i)/(aMotorPos.getDataSize() * aSlList.getDataSize())));
-        }
-    }
-
-    thread.join();
-    return Env;
+    // thread.join();
+    RECON_INFO("Reflect reconstructing done");
+    return Matrix();
 }
--- a/src/startReconstructions.cpp
+++ b/src/startReconstructions.cpp
@@ -4,13 +4,13 @@
 #include "MatlabWriter.h"
 #include "config/config.h"
 #include "log/log.h"
-#include "common/DICOMExporter.h"
 #include "common/getMeasurementMetaData.h"
 #include "common/getGeometryInfo.h"
 #include "common/estimatePulseLength.h"
 #include "log/notify.h"
 #include "transmissionReconstruction/startTransmissionReconstruction.h"

+#include <omp.h>
 #include <string>
 #include <iostream>
 #include "Parser/Parser.h"
@@ -33,11 +33,12 @@ using namespace Aurora;

 int Recon::startReconstructions(const std::string& aDataPath, const std::string& aDataRefPath, const std::string& aOutputPath)
 {
-    MatlabWriter writer(aOutputPath);
+    // omp_set_num_threads(32);
+    // mkl_set_num_threads(64);
+
+    // MatlabWriter writer(aOutputPath);
    Parser dataParser(aDataPath);
    Parser refParser(aDataRefPath);
-    Recon::DICOMExporter exporter(dataParser.getPatientData());
-    exporter.setExportBasePath(aOutputPath);
    if(!dataParser.getShotList()->isValid())
    {
        RECON_INFO("startReconstructions with {0}.",aDataPath);
@@ -224,8 +225,6 @@ int Recon::startReconstructions(const std::string& aDataPath, const std::string&
        transmissionResult = startTransmissionReconstruction(mp_inter, mpRef_inter, slList_inter, snList_inter, rlList_inter, rnList_inter, temp, tempRef, geom, geomRef, expInfo, expInfoRef, preComputes, &dataParser, &refParser);
        attAvailable = true;
        sosAvailable = true;
-        exporter.exportDICOM(transmissionResult.recoSOS, DICOMExporter::SOS);
-        exporter.exportDICOM(transmissionResult.recoATT, DICOMExporter::ATT);
        //  writer.setMatrix(transmissionResult.recoSOS, "sos");
        //  writer.setMatrix(transmissionResult.recoATT, "att");
    }
@@ -255,8 +254,7 @@ int Recon::startReconstructions(const std::string& aDataPath, const std::string&
        //Recon::notifyProgress(25);
        RECON_INFO("Start reflectionRecostruction.");
        Matrix env = startReflectionReconstruction(&dataParser, preProcessData.saftMode, mp_inter, slList_inter, snList_inter, rlList_inter, rnList_inter, geom, preProcessData.transRecos, expInfo, preComputes);
-         writer.setMatrix(env, "reflect");
-        exporter.exportDICOM(env, Recon::DICOMExporter::REFL);
+        //  writer.setMatrix(env, "reflect");
        // Recon::notifyProgress(99);
    }
        //  writer.write();
--- a/src/transmissionReconstruction/detection/detection.cpp
+++ b/src/transmissionReconstruction/detection/detection.cpp
@@ -12,8 +12,8 @@

 #include "common/getMeasurementMetaData.h"
 #include "config/config.h"
-#include "calculateBankDetectAndHilbertTransformation.hpp"
-#include "transmissionReconstruction/detection/detection.cuh"
+
+#include "log/log.h"

 using namespace Aurora;
 namespace Recon {
@@ -156,6 +156,7 @@ namespace Recon {
            startPosRef[i] =  std::floor(std::max( tof2[i],(float)1.0));
            endPosRef[i] = std::ceil(std::min(sizeAscan[1], tof2[i]+detectionWindowATT));
        }
+        RECON_INFO("calculateAttenuation start");
        return calculateAttenuation(envelopeOfAScan,startPos,endPos,envelopeOfReferenceAScan,startPosRef,endPosRef);
    }

@@ -259,6 +260,7 @@ namespace Recon {
        TimeWindowResult timeResult2;
        timeResult2.AscanBlockProcessed = AscanRefBlock;
        if (useTimeWindowing == 1) {
+            RECON_INFO("applyTimeWindowing start");
            timeResult1 = applyTimeWindowing(
                AscanBlock, sampleRate, distBlock, sosWaterBlock,
                expectedSOSWater, offsetElectronicSamples, detectionWindowSOS,
@@ -269,6 +271,7 @@ namespace Recon {
                minSpeedOfSound, maxSpeedOfSound, gaussWindow);

            diffStartSearch = timeResult1.startSearch - timeResult2.startSearch;
+            RECON_INFO("applyTimeWindowing end");
        }
        auto _AscanBlock = timeResult1.AscanBlockProcessed;
        auto _AscanRefBlock = timeResult2.AscanBlockProcessed;
@@ -279,12 +282,12 @@ namespace Recon {
        auto mxl = std::min(maxlag, m - 1);
        auto ceilLog2 = nextpow2(2 * m - 1);
        auto m2 = pow(2, ceilLog2);
-
+        RECON_INFO("fft start");
        auto x = fft(_AscanBlock, m2);
        auto y = fft(_AscanRefBlock, m2);
        auto c_1_1 = x * conj(y);
        auto c_1_2 = ifft(c_1_1);
-
+        RECON_INFO("fft end");
        auto c1 = real(c_1_2);
        auto c = zeros(mxl + mxl + 1, c1.getDimSize(1));
        c.setBlock(0, 0, mxl-1, c1.block(0, m2-mxl,m2-1));
@@ -311,11 +314,13 @@ namespace Recon {
        DetectResult result;
        result.sosValue = sosValue;
        auto tofRel = tof - distBlock / sosWaterBlock;
+        RECON_INFO("detectAttVectorized start");
        result.att = detectAttVectorized(
            _AscanBlock, _AscanRefBlock, distRefBlock, sosWaterRefBlock,
            tof, aScanReconstructionFrequency, offsetElectronic,
            detectionWindowATT);
        result.tof = tofRel;
+        RECON_INFO("detectAttVectorized end");
        return result;
    }

@@ -367,9 +372,6 @@ namespace Recon {
        resDetect = Aurora::zeros(1,N);
        resEnvelopeD = Aurora::malloc(totalSize);
        resEnvelopeRefD = Aurora::malloc(totalSize);
-        calculateBankDetectAndHilbertTransformation(
-                _AscanBlock_trim.getData(), _AscanRefBlock_trim.getData(), N, M, resampleFactor, nthreads,
-                resDetect.getData(), resEnvelopeD, resEnvelopeRefD);
        auto resEnvelope =Matrix::New(resEnvelopeD,M,N);
        auto resEnvelopeRef =Matrix::New(resEnvelopeRefD,M,N);
        //floor(size(AscanBlock, 1)*inits.resampleFactor / 2 - 1),
@@ -415,15 +417,15 @@ namespace Recon {
        return result;
    }

-    DetectResult transmissionDetection(const Aurora::CudaMatrix &AscanBlock,
-                                       const Aurora::CudaMatrix &AscanRefBlock,
-                                       const Aurora::CudaMatrix &distBlock,
-                                       const Aurora::CudaMatrix &distRefBlock,
+    DetectResult transmissionDetection(const Aurora::Matrix &AscanBlock,
+                                       const Aurora::Matrix &AscanRefBlock,
+                                       const Aurora::Matrix &distBlock,
+                                       const Aurora::Matrix &distRefBlock,
                                       const Aurora::Matrix &sosWaterBlock,
                                       const Aurora::Matrix &sosWaterRefBlock,
                                       float expectedSOSWater) {
-        auto _sosWaterBlock = temperatureToSoundSpeed(sosWaterBlock, "marczak").toDeviceMatrix();
-        auto _sosWaterRefBlock = temperatureToSoundSpeed(sosWaterRefBlock, "marczak").toDeviceMatrix();
+        auto _sosWaterBlock = temperatureToSoundSpeed(sosWaterBlock, "marczak");
+        auto _sosWaterRefBlock = temperatureToSoundSpeed(sosWaterRefBlock, "marczak");
        switch (Recon::transParams::version) {
        // case 1: {
        //     return detectTofAndAttMex(
@@ -444,9 +446,9 @@ namespace Recon {
                Recon::transParams::offsetElectronic, Recon::transParams::detectionWindowSOS, Recon::transParams::minSpeedOfSound,
                Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow);
            DetectResult ret;
-            ret.att = r.att.toHostMatrix();
-            ret.sosValue = r.sosValue.toHostMatrix();
-            ret.tof = r.tof.toHostMatrix();
+            ret.att = r.att;
+            ret.sosValue = r.sosValue;
+            ret.tof = r.tof;
            return ret;
        }
    }
--- a/src/transmissionReconstruction/detection/detection.cuh_
+++ b/src/transmissionReconstruction/detection/detection.cuh_
--- a/src/transmissionReconstruction/detection/detection.h
+++ b/src/transmissionReconstruction/detection/detection.h
@@ -75,8 +75,8 @@ detectTofAndAttMex(

 DetectResult 
 transmissionDetection(
-    const Aurora::CudaMatrix &AscanBlock, const Aurora::CudaMatrix &AscanRefBlock,
-    const Aurora::CudaMatrix &distBlock, const Aurora::CudaMatrix &distRefBlock,
+    const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
+    const Aurora::Matrix &distBlock, const Aurora::Matrix &distRefBlock,
    const Aurora::Matrix &sosWaterBlock, const Aurora::Matrix &sosWaterRefBlock, float expectedSOSWater);

 } // namespace Recon
--- a/src/transmissionReconstruction/detection/getTransmissionData.cpp
+++ b/src/transmissionReconstruction/detection/getTransmissionData.cpp
@@ -1,11 +1,9 @@
 #include "getTransmissionData.h"
-#include "getTransmissionData.cuh"
-#include "AuroraDefs.h"
-#include "CudaMatrix.h"
 #include "Function.h"
 #include "Function1D.h"
 #include "Function2D.h"
 #include "common/dataBlockCreation/removeDataFromArrays.h"
+#include "log/log.h"
 #include "log/notify.h"
 #include "src//config/config.h"
 #include "src/common/getGeometryInfo.h"
@@ -31,11 +29,6 @@
 #include <mutex>
 #include <condition_variable>

-#include <cuda_runtime.h>
-#include "Function1D.cuh"
-#include "Function2D.cuh"
-#include <sys/time.h>
-
 using namespace Recon;
 using namespace Aurora;

@@ -50,8 +43,8 @@ namespace
        Matrix waterTempBlock;
        MetaInfos metaInfos;

-        CudaMatrix ascanBlock;
-        CudaMatrix ascanBlockRef;
+        Matrix ascanBlock;
+        Matrix ascanBlockRef;
        Matrix dists;
        Matrix distRefBlock;
        Matrix waterTempRefBlock;
@@ -64,59 +57,35 @@ namespace
    std::mutex PROCESS_BUFFER_MUTEX;
    std::condition_variable PROCESS_BUFFER_CONDITION;
    int BUFFER_COUNT = 0;
-    int BUFFER_SIZE = 4;//<=8
+    int BUFFER_SIZE = 1;

-void printTime()
-{
-    struct timeval tpend;
-    gettimeofday(&tpend,NULL);
-    int secofday = (tpend.tv_sec + 3600 * 8 ) % 86400;
-    int hours = secofday / 3600;
-    int minutes = (secofday - hours * 3600 ) / 60;
-    int seconds = secofday % 60; 
-    int milliseconds = tpend.tv_usec/1000;
-    std::cout<< hours << ":" <<minutes<<":"<<seconds<<"."<<milliseconds<<std::endl;
-}
-
-    CudaMatrix prepareAScansForTransmissionDetection(const CudaMatrix& aAscanBlock, const CudaMatrix& aGainBlock)
+    Matrix prepareAScansForTransmissionDetection(const Matrix& aAscanBlock, const Matrix& aGainBlock)
    {
-        CudaMatrix result = aAscanBlock / repmat(aGainBlock, aAscanBlock.getDimSize(0), 1);
+        RECON_INFO("prepareAScansForTransmissionDetection start");
+        Matrix result = aAscanBlock / repmat(aGainBlock, aAscanBlock.getDimSize(0), 1);
        result = result - repmat(mean(result,FunctionDirection::Column), result.getDimSize(0), 1);
+        RECON_INFO("prepareAScansForTransmissionDetection end");
        return result;
    }

-    Aurora::CudaMatrix calculateSnr(const Aurora::CudaMatrix &aMDataBlock,
-                                float aReferenceNoise) {
-        auto maxSignal = max(abs(aMDataBlock));
-        auto snrBlock = 10 * log(maxSignal / aReferenceNoise, 10);
-        return snrBlock;
-    }
-
    BlockOfTransmissionData getBlockOfTransmissionData(const Matrix& aMp, const Matrix& aMpRef, const Matrix& aSl, const Matrix& aSn, const Matrix& aRlList, const Matrix& aRnList,
                                    const TasTemps& aTasTemps, const Matrix& aExpectedSOSWater, GeometryInfo aGeom, GeometryInfo& aGeomRef,
                                    const Matrix& aSnrRmsNoise, const Matrix& aSnrRmsNoiseRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
                                    const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
    {
        BlockOfTransmissionData result;
-        MetaInfos metaInfos;        
-//printTime();
-//2500ms        
-        auto blockData = getAscanBlockPreprocessedCuda(aParser, aMp, aSl, aSn, aRlList, aRnList, aGeom, aExpInfo, true, true);
-        auto blockDataRef = getAscanBlockPreprocessedCuda(aParserRef, aMpRef, aSl, aSn, aRlList, aRnList, aGeomRef, aExpInfoRef, true, true);
-//printTime();
-//180ms       
-    //    auto t1 = blockData.ascanBlockPreprocessed.toDeviceMatrix();
-    //    auto t2 = blockDataRef.ascanBlockPreprocessed.toDeviceMatrix();
-    //     auto t3 = blockData.gainBlock.toDeviceMatrix();
-    //    auto t4 = blockDataRef.gainBlock.toDeviceMatrix();
-//printTime();       
-//20ms       
-        CudaMatrix ascanBlock = prepareAScansForTransmissionDetection(blockData.ascanBlockPreprocessed,blockData.gainBlock);
-        CudaMatrix ascanBlockRef = prepareAScansForTransmissionDetection(blockDataRef.ascanBlockPreprocessed,blockDataRef.gainBlock);
-//printTime();
-//20ms
-        blockData.ascanBlockPreprocessed = CudaMatrix();
-        blockDataRef.ascanBlockPreprocessed = CudaMatrix();
+        MetaInfos metaInfos;
+        RECON_INFO("getAscanBlockPreprocessed start");
+        auto blockData = getAscanBlockPreprocessed(aParser, aMp, aSl, aSn, aRlList, aRnList, aGeom, aExpInfo, true, true);
+        auto blockDataRef = getAscanBlockPreprocessed(aParserRef, aMpRef, aSl, aSn, aRlList, aRnList, aGeomRef, aExpInfoRef, true, true);
+        RECON_INFO("getAscanBlockPreprocessed end");
+        RECON_INFO("prepareAScansForTransmissionDetection start");
+        Matrix ascanBlock = prepareAScansForTransmissionDetection(blockData.ascanBlockPreprocessed, blockData.gainBlock);
+        Matrix ascanBlockRef = prepareAScansForTransmissionDetection(blockDataRef.ascanBlockPreprocessed, blockDataRef.gainBlock);
+        RECON_INFO("prepareAScansForTransmissionDetection end");
+        blockData.ascanBlockPreprocessed = Matrix();
+        blockDataRef.ascanBlockPreprocessed = Matrix();
+        RECON_INFO("Matrix * + start");
        if(aExpInfo.Hardware == "USCT3dv3")
        {
            Matrix channelList = precalculateChannelList(aRlList, aRnList, aExpInfo, aPreComputes);
@@ -132,48 +101,72 @@ void printTime()
                channelListBlockData[i] = channelList[ind[i] - 1];
            }
            Matrix channelListBlock = Matrix::New(channelListBlockData, 1, channelListBlockSize);
-//printTime();  
-//20ms          
-            CudaMatrix fx = fft(ascanBlock);
-//printTime();
-//50ms 
-            // float* fhData =  nullptr;
-            // cudaMalloc((void**)&fhData, sizeof(float) * aExpInfo.matchedFilter.getDimSize(0) * channelListBlockSize * Aurora::Complex);
-            // CudaMatrix fh = CudaMatrix::fromRawData(fhData, aExpInfo.matchedFilter.getDimSize(0), channelListBlockSize, 1, Aurora::Complex);
-            // size_t matchedFilterRowDataSize = aExpInfo.matchedFilter.getDimSize(0)*2;
-            // for(size_t i=0; i<channelListBlockSize; ++i)
-            // {
-            //     cudaMemcpy(fhData, aExpInfo.matchedFilter.getData() + (size_t)(channelListBlock[i] - 1) * matchedFilterRowDataSize, sizeof(float) * matchedFilterRowDataSize, cudaMemcpyHostToDevice);
-            //     fhData += matchedFilterRowDataSize;
-            // }
-            Aurora::CudaMatrix matchedFilterDevice = aExpInfo.matchedFilter.toDeviceMatrix();
-            Aurora::CudaMatrix channelListBlockDevice = channelListBlock.toDeviceMatrix();
-            Aurora::CudaMatrix fh = createFhMatrix(matchedFilterDevice, channelListBlockDevice);        
- //printTime();
- //20ms
-            CudaMatrix complex = getTransmissionDataSubFunction(fx, fh);
+            Matrix fx = fft(ascanBlock);
+            float* fhData =  Aurora::malloc(aExpInfo.matchedFilter.getDimSize(0) * channelListBlockSize, true);
+            Matrix fh = Matrix::New(fhData, aExpInfo.matchedFilter.getDimSize(0), channelListBlockSize, 1, Aurora::Complex);
+            size_t matchedFilterRowDataSize = aExpInfo.matchedFilter.getDimSize(0)*2;
+            for(size_t i=0; i<channelListBlockSize; ++i)
+            {
+                cblas_scopy(matchedFilterRowDataSize, aExpInfo.matchedFilter.getData() + (size_t)(channelListBlock[i] - 1) * matchedFilterRowDataSize, 1 , fhData ,1);
+                fhData += matchedFilterRowDataSize;
+            }
+            // Matrix fxReal = Aurora::real(fx);
+            // Matrix fhReal = Aurora::real(fh);
+            // Matrix fxImag = Aurora::imag(fx);
+            // Matrix fhImag = Aurora::imag(fh);
+            // Matrix real = fxReal * fhReal + fxImag * fhImag;
+            // Matrix image =  fxImag * fhReal - fxReal * fhImag;
+            float* value1 = Aurora::malloc(fx.getDataSize());
+            vsMulI(fx.getDataSize(), fx.getData(), 2, fh.getData(), 2, value1, 1);
+            float* value2 = Aurora::malloc(fx.getDataSize());
+            vsMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData() + 1, 2, value2, 1);
+            float* realData = Aurora::malloc(fx.getDataSize());
+            vsAdd(fx.getDataSize(), value1, value2, realData);
+            Matrix real = Matrix::New(realData, fx.getDimSize(0), fx.getDimSize(1));
+
+            vsMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData(), 2, value1, 1);
+            vsMulI(fx.getDataSize(), fx.getData(), 2, fh.getData() + 1, 2, value2, 1);
+            float* imagData = Aurora::malloc(fx.getDataSize());
+            vsSub(fx.getDataSize(), value1, value2, imagData);
+            Matrix image = Matrix::New(imagData, fx.getDimSize(0), fx.getDimSize(1));
+
+            float* complexData = Aurora::malloc(real.getDataSize(), true);
+            cblas_scopy(real.getDataSize(), real.getData(), 1 , complexData ,2);
+            cblas_scopy(image.getDataSize(), image.getData(), 1 , complexData + 1 ,2);
+            Matrix complex = Matrix::New(complexData, real.getDimSize(0), real.getDimSize(1), 1, Aurora::Complex);
            ascanBlock = Aurora::real(ifft(complex));
-//printTime();            
-//20s
+
            fx = fft(ascanBlockRef);
- //printTime(); 
-//50ms            
-            // cudaMalloc((void**)&fhData, sizeof(float) * aExpInfo.matchedFilter.getDimSize(0) * channelListBlockSize * Aurora::Complex);
-            // fh = CudaMatrix::fromRawData(fhData, aExpInfoRef.matchedFilter.getDimSize(0), channelListBlockSize, 1, Aurora::Complex);
-            // matchedFilterRowDataSize = aExpInfoRef.matchedFilter.getDimSize(0)*2;
-            // for(size_t i=0; i<channelListBlockSize; ++i)
-            // {
-            //     cudaMemcpy(fhData, aExpInfoRef.matchedFilter.getData() + (size_t)(channelListBlock[i] - 1) * matchedFilterRowDataSize, sizeof(float) * matchedFilterRowDataSize, cudaMemcpyHostToDevice);
-            //     fhData += matchedFilterRowDataSize;
-            // }
-            matchedFilterDevice = aExpInfoRef.matchedFilter.toDeviceMatrix();
-            fh = createFhMatrix(matchedFilterDevice, channelListBlockDevice);  
- //printTime();
- //20ms
-            complex = getTransmissionDataSubFunction(fx, fh);
+            fhData =  Aurora::malloc(aExpInfoRef.matchedFilter.getDimSize(0) * channelListBlockSize, true);
+            fh = Matrix::New(fhData, aExpInfoRef.matchedFilter.getDimSize(0), channelListBlockSize, 1, Aurora::Complex);
+            matchedFilterRowDataSize = aExpInfoRef.matchedFilter.getDimSize(0)*2;
+            for(size_t i=0; i<channelListBlockSize; ++i)
+            {
+                cblas_scopy(matchedFilterRowDataSize, aExpInfoRef.matchedFilter.getData() + (size_t)(channelListBlock[i] - 1) * matchedFilterRowDataSize, 1 , fhData ,1);
+                fhData += matchedFilterRowDataSize;
+            }
+            //  real = Aurora::real(fx) * Aurora::real(fh) + Aurora::imag(fx) * Aurora::imag(fh);
+            //  image =  Aurora::imag(fx) * Aurora::real(fh) - Aurora::real(fx) * Aurora::imag(fh);
+            vsMulI(fx.getDataSize(), fx.getData(), 2, fh.getData(), 2, value1, 1);
+            vsMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData() + 1, 2, value2, 1);
+            realData = Aurora::malloc(fx.getDataSize());
+            vsAdd(fx.getDataSize(), value1, value2, realData);
+            real = Matrix::New(realData, fx.getDimSize(0), fx.getDimSize(1));
+
+            vsMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData(), 2, value1, 1);
+            vsMulI(fx.getDataSize(), fx.getData(), 2, fh.getData() + 1, 2, value2, 1);
+            imagData = Aurora::malloc(fx.getDataSize());
+            vsSub(fx.getDataSize(), value1, value2, imagData);
+            image = Matrix::New(imagData, fx.getDimSize(0), fx.getDimSize(1));
+            Aurora::free(value1);
+            Aurora::free(value2);
+
+            complexData = Aurora::malloc(real.getDataSize(), true);
+            cblas_scopy(real.getDataSize(), real.getData(), 1 , complexData ,2);
+            cblas_scopy(image.getDataSize(), image.getData(), 1 , complexData + 1 ,2);
+            complex = Matrix::New(complexData, real.getDimSize(0), real.getDimSize(1), 1, Aurora::Complex);
            ascanBlockRef = Aurora::real(ifft(complex));
-//printTime();
-//20ms
+            RECON_INFO("Matrix * + end");
        }
        else
        {
@@ -183,27 +176,23 @@ void printTime()

        if(transParams::applyCalib)
        {
-            metaInfos.snrValues = calculateSnr(ascanBlock, aSnrRmsNoise[0]).toHostMatrix();
-            metaInfos.snrValuesRef = calculateSnr(ascanBlockRef, aSnrRmsNoiseRef[0]).toHostMatrix();
+            RECON_INFO("calculateSnr start");
+            metaInfos.snrValues = calculateSnr(ascanBlock, aSnrRmsNoise[0]);
+            metaInfos.snrValuesRef = calculateSnr(ascanBlockRef, aSnrRmsNoiseRef[0]);
+            RECON_INFO("calculateSnr end");
        }
-// printTime();
- //3ms
        Matrix dists = distanceBetweenTwoPoints(blockData.senderPositionBlock, blockData.receiverPositionBlock);
        Matrix distRefBlock = distanceBetweenTwoPoints(blockDataRef.senderPositionBlock, blockDataRef.receiverPositionBlock);
- //printTime();
- //2ms
        Matrix waterTempBlock = calculateWaterTemperature(aTasTemps.waterTempPreCalc_sl, aTasTemps.waterTempPreCalc_rl, blockData.slBlock, blockData.rlBlock, blockData.mpBlock);
        Matrix waterTempRefBlock  = calculateWaterTemperature(aTasTemps.waterTempRefPreCalc_sl, aTasTemps.waterTempRefPreCalc_rl, blockData.slBlock, blockData.rlBlock, blockDataRef.mpBlock);
- // printTime();
- // 1ms       
-        // if(transParams::saveDetection || transParams::outlierOnTasDetection || transParams::saveDebugInfomation)
-        // {
-        //     metaInfos.mpBlock = blockData.mpBlock;
-        //     metaInfos.slBlock = blockData.slBlock;
-        //     metaInfos.snBlock = blockData.snBlock;
-        //     metaInfos.rlBlock = blockData.rlBlock;
-        //     metaInfos.rnBlock = blockData.rnBlock;
-        // }
+        if(transParams::saveDetection || transParams::outlierOnTasDetection || transParams::saveDebugInfomation)
+        {
+            metaInfos.mpBlock = blockData.mpBlock;
+            metaInfos.slBlock = blockData.slBlock;
+            metaInfos.snBlock = blockData.snBlock;
+            metaInfos.rlBlock = blockData.rlBlock;
+            metaInfos.rnBlock = blockData.rnBlock;
+        }
        result.metaInfos = metaInfos;
        result.senderBlock = blockData.senderPositionBlock;
        result.receiverBlock = blockData.receiverPositionBlock;
@@ -218,7 +207,7 @@ void printTime()
        // DetectResult detect = transmissionDetection(ascanBlock, ascanBlockRef, dists, distRefBlock, waterTempBlock, waterTempRefBlock, aExpectedSOSWater[0]);
        // result.attData = detect.att;
        // result.tofData = detect.tof;
-//printTime();
+        
        return result;
    }

@@ -227,17 +216,18 @@ void printTime()
 void  getBlockOfTransmissionDataInThread(size_t aIndex, const Matrix& aMp, const Matrix& aMpRef, const Matrix& aSl, const Matrix& aSn, const Matrix& aRlList, const Matrix& aRnList,
                                         const TasTemps& aTasTemps, const Matrix& aExpectedSOSWater, GeometryInfo aGeom, GeometryInfo aGeomRef,
                                         const Matrix& aSnrRmsNoise, const Matrix& aSnrRmsNoiseRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
-                                         const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef, unsigned int aGPUId)
+                                         const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
 {
-    cudaSetDevice(aGPUId); 
    auto buffer = getBlockOfTransmissionData(aMp, aMpRef, aSl, aSn, aRlList, aRnList, aTasTemps, 
                                             aExpectedSOSWater, aGeom, aGeomRef, aSnrRmsNoise, aSnrRmsNoiseRef,
                                             aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef);
    std::unique_lock<std::mutex> lock(CREATE_BUFFER_MUTEX); 
    BLOCK_OF_TRANSIMISSIONDARA_BUFFER[std::to_string(aIndex)] = buffer;
-    std::cout<<"Add: "<<aIndex<<std::endl;
+    RECON_INFO("add");
+    std::cout<<"add: "<<aIndex<<std::endl;
    lock.unlock();
    PROCESS_BUFFER_CONDITION.notify_one();
+    
 }

 void  createThreadForGetBlockOfTransmissionData(const Matrix& aMotorPos, const Matrix& aMotoPosRef, const Matrix& aSlList, const Matrix& aSnList, const Matrix& aRlList, const Matrix& aRnList,
@@ -246,14 +236,20 @@ void  createThreadForGetBlockOfTransmissionData(const Matrix& aMotorPos, const M
                                         const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
 {
    size_t vectorSize = aMotorPos.getDataSize() * (aSlList.getDataSize() / transParams::senderTASSize) * (aSnList.getDataSize() / transParams::senderElementSize);
-    std::thread speedUpThread[vectorSize];
+    std::thread speedUpThread[/*vectorSize*/1];

    for(int i=0; i<aMotorPos.getDataSize(); ++i)
    {
-        for(int j=0; j<aSlList.getDataSize() / transParams::senderTASSize; ++j)
+        for(int j=0; j</*aSlList.getDataSize() / transParams::senderTASSize*/1; ++j)
        {
-            for(int k=0; k<aSnList.getDataSize() / transParams::senderElementSize; ++k)
+            for(int k=0; k</*aSnList.getDataSize() / transParams::senderElementSize*/1; ++k)
            {
+    // for(int i=0; i<1; ++i)
+    // {
+    //     for(int j=0; j<1; ++j)
+    //     {
+    //         for(int k=0; k<3; ++k)
+    //         {
                size_t index = i * (aSlList.getDataSize() / transParams::senderTASSize) * (aSnList.getDataSize() / transParams::senderElementSize) + 
                               j * (aSnList.getDataSize() / transParams::senderElementSize) + k;
                Matrix mp = aMotorPos(i).toMatrix();
@@ -264,8 +260,8 @@ void  createThreadForGetBlockOfTransmissionData(const Matrix& aMotorPos, const M
                std::unique_lock<std::mutex> lock(CREATE_BUFFER_MUTEX); 
                CREATE_BUFFER_CONDITION.wait(lock, []{return BUFFER_COUNT<BUFFER_SIZE;});
                ++BUFFER_COUNT;
-                lock.unlock();
-                speedUpThread[index] =  std::thread(getBlockOfTransmissionDataInThread,index,mp,mpRef,sl,sn,aRlList,aRnList,aTasTemps,aExpectedSOSWater,aGeom,aGeomRef,aSnrRmsNoise,aSnrRmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef, index % BUFFER_SIZE);
+                lock.unlock(); 
+                speedUpThread[index] =  std::thread(getBlockOfTransmissionDataInThread,index,mp,mpRef,sl,sn,aRlList,aRnList,aTasTemps,aExpectedSOSWater,aGeom,aGeomRef,aSnrRmsNoise,aSnrRmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef);
            }
        }
    }
@@ -289,12 +285,12 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
    //     % Load transmission detection data
    //     writeReconstructionLog('Loading transmission detection data. All available data from given motor positions are taken.', 1);
    // [tofDataTotal, attDataTotal, senderList, receiverList, waterTempList, dataInfo] = loadTransmissionDetectionData(transParams.pathSaveDetection, transParams.pathData, motorPos, expInfo.rootMeasUniqueID);
+    
    TasTemps tasTemps;
    tasTemps.waterTempPreCalc_rl = extractTasTemperature(aTemp.tasTemperature, aRlList, aMotorPos, aTemp.jumoTemp, transParams::minTemperature, transParams::maxTemperature);
    tasTemps.waterTempPreCalc_sl = extractTasTemperature(aTemp.tasTemperature, aSlList, aMotorPos, aTemp.jumoTemp, transParams::minTemperature, transParams::maxTemperature);
    tasTemps.waterTempRefPreCalc_rl = extractTasTemperature(aTempRef.tasTemperature, aRlList, aMotoPosRef, aTempRef.jumoTemp, transParams::minTemperature, transParams::maxTemperature);
    tasTemps.waterTempRefPreCalc_sl = extractTasTemperature(aTempRef.tasTemperature, aSlList, aMotoPosRef, aTempRef.jumoTemp, transParams::minTemperature, transParams::maxTemperature);
-
    aGeom.sensData = precalcSensitivity(aSlList, aSnList, aRlList, aRnList, aMotorPos, aGeom);
    aGeomRef.sensData = aGeom.sensData;

@@ -304,7 +300,6 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
        rmsNoise = estimateNoiseValueFromAScans(aSnList, aRnList, aGeom, aExpInfo, aParser);
        rmsNoiseRef = estimateNoiseValueFromAScans(aSnList, aRnList, aGeom, aExpInfo, aParserRef);
    }
-
    size_t numScans = aMotorPos.getDataSize() * aSlList.getDataSize() * aSnList.getDataSize() * aRlList.getDataSize() * aRnList.getDataSize();
    Matrix tofDataTotal = Matrix::fromRawData(new float[numScans], 1, numScans) + NAN;
    Matrix attDataTotal = Matrix::fromRawData(new float[numScans], 1, numScans) + NAN;
@@ -327,7 +322,6 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
        rlBlockTotal = zeros(1,numScans,1);
        rnBlockTotal = zeros(1,numScans,1);
    }
-
    std::thread speedUpThread = std::thread(createThreadForGetBlockOfTransmissionData,aMotorPos,aMotoPosRef,aSlList,aSnList,aRlList,aRnList,tasTemps,aTemp.expectedSOSWater,aGeom,aGeomRef,rmsNoise,rmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef);
    int numData = 0;
    int numPossibleScans = 0;
@@ -339,13 +333,19 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
    {
        std::cerr << "Failed to set thread priority" << std::endl;
    }
-
+    int  finish_count=0;
    for(int i=0; i<aMotorPos.getDataSize(); ++i)
    {
-        for(int j=0; j<aSlList.getDataSize() / transParams::senderTASSize; ++j)
+        for(int j=0; j</*aSlList.getDataSize() / transParams::senderTASSize*/1; ++j)
        {
-            for(int k=0; k<aSnList.getDataSize() / transParams::senderElementSize; ++k)
+            for(int k=0; k</*aSnList.getDataSize() / transParams::senderElementSize*/1; ++k)
            {
+    // for(int i=0; i<1; ++i)
+    // {
+    //     for(int j=0; j<1; ++j)
+    //     {
+    //         for(int k=0; k<3; ++k)
+    //         {
                size_t index = i * (aSlList.getDataSize() / transParams::senderTASSize) * (aSnList.getDataSize() / transParams::senderElementSize) + 
                                j * (aSnList.getDataSize() / transParams::senderElementSize) + k;
                std::unique_lock<std::mutex> lock(PROCESS_BUFFER_MUTEX);
@@ -354,15 +354,17 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
                lock.unlock();                

                auto blockData = BLOCK_OF_TRANSIMISSIONDARA_BUFFER[std::to_string(index)];
-                cudaSetDevice(index % BUFFER_SIZE); 
+                RECON_INFO("transmissionDetection start");
                DetectResult detect = transmissionDetection( blockData.ascanBlock, blockData.ascanBlockRef, 
-                                                             blockData.dists.toDeviceMatrix(), blockData.distRefBlock.toDeviceMatrix(), 
+                                                             blockData.dists, blockData.distRefBlock, 
                                                             blockData.waterTempBlock, blockData.waterTempRefBlock,
                                                             aTemp.expectedSOSWater[0]);
+                RECON_INFO("transmissionDetection end");                                                             
                blockData.attData = detect.att;
                blockData.tofData = detect.tof;
                BlockOfTransmissionData transmissionBlock=blockData;
                size_t numUsedData = transmissionBlock.senderBlock.getDimSize(1);
+                RECON_INFO("cblas_scopy start");
                if(transParams::applyCalib)
                {
                    cblas_scopy(numUsedData, transmissionBlock.metaInfos.snrValues.getData(), 1, snrValues.getData() + numData, 1);
@@ -376,24 +378,29 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
                cblas_scopy(numUsedData, transmissionBlock.attData.getData(), 1, attDataTotal.getData() + numData, 1);
                cblas_scopy(numUsedData, transmissionBlock.waterTempBlock.getData(), 1, waterTempList.getData() + numData, 1);
                
-                // if(transParams::saveDetection || transParams::outlierOnTasDetection || transParams::saveDebugInfomation)
-                // {
-                //     cblas_scopy(numUsedData, transmissionBlock.metaInfos.mpBlock.getData(), 1, mpBlockTotal.getData() + numData, 1);
-                //     cblas_scopy(numUsedData, transmissionBlock.metaInfos.slBlock.getData(), 1, slBlockTotal.getData() + numData, 1);
-                //     cblas_scopy(numUsedData, transmissionBlock.metaInfos.snBlock.getData(), 1, snBlockTotal.getData() + numData, 1);
-                //     cblas_scopy(numUsedData, transmissionBlock.metaInfos.rlBlock.getData(), 1, rlBlockTotal.getData() + numData, 1);
-                //     cblas_scopy(numUsedData, transmissionBlock.metaInfos.rnBlock.getData(), 1, rnBlockTotal.getData() + numData, 1);
-                // }
+                if(transParams::saveDetection || transParams::outlierOnTasDetection || transParams::saveDebugInfomation)
+                {
+                    cblas_scopy(numUsedData, transmissionBlock.metaInfos.mpBlock.getData(), 1, mpBlockTotal.getData() + numData, 1);
+                    cblas_scopy(numUsedData, transmissionBlock.metaInfos.slBlock.getData(), 1, slBlockTotal.getData() + numData, 1);
+                    cblas_scopy(numUsedData, transmissionBlock.metaInfos.snBlock.getData(), 1, snBlockTotal.getData() + numData, 1);
+                    cblas_scopy(numUsedData, transmissionBlock.metaInfos.rlBlock.getData(), 1, rlBlockTotal.getData() + numData, 1);
+                    cblas_scopy(numUsedData, transmissionBlock.metaInfos.rnBlock.getData(), 1, rnBlockTotal.getData() + numData, 1);
+                }
+                RECON_INFO("cblas_scopy end");
                numData += numUsedData;
                std::unique_lock<std::mutex> lockBufferCount(CREATE_BUFFER_MUTEX); 
                BLOCK_OF_TRANSIMISSIONDARA_BUFFER.erase(std::to_string(index));
                --BUFFER_COUNT;
+                finish_count++;
                lockBufferCount.unlock();
+                RECON_INFO("Remove");
                std::cout<<"Remove: "<<index<<std::endl;
                CREATE_BUFFER_CONDITION.notify_one();
            }
-            //Recon::notifyProgress(6+10*((i+1)*(j+1)/(aMotorPos.getDataSize()*(aSlList.getDataSize()/ transParams::senderTASSize))));   
+            Recon::notifyProgress(6+10*((i+1)*(j+1)/(aMotorPos.getDataSize()*(aSlList.getDataSize()/ transParams::senderTASSize))));   
        }
+        std::cout<<"Transmission reconstruct done! "<< std::endl;
+        std::cout.flush();
    }
    speedUpThread.join();

@@ -420,14 +427,14 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
    tofDataTotal = removeDataFromArrays(tofDataTotal, filter);
    attDataTotal = removeDataFromArrays(attDataTotal, filter);
    waterTempList = removeDataFromArrays(waterTempList, filter);
-    // if(transParams::saveDebugInfomation || transParams::outlierOnTasDetection || transParams::saveDetection)
-    // {
-    //     mpBlockTotal = removeDataFromArrays(mpBlockTotal, filter);
-    //     slBlockTotal = removeDataFromArrays(slBlockTotal, filter);
-    //     snBlockTotal = removeDataFromArrays(snBlockTotal, filter);
-    //     rlBlockTotal = removeDataFromArrays(rlBlockTotal, filter);
-    //     rnBlockTotal = removeDataFromArrays(rnBlockTotal, filter);
-    // }
+    if(transParams::saveDebugInfomation || transParams::outlierOnTasDetection || transParams::saveDetection)
+    {
+        mpBlockTotal = removeDataFromArrays(mpBlockTotal, filter);
+        slBlockTotal = removeDataFromArrays(slBlockTotal, filter);
+        snBlockTotal = removeDataFromArrays(snBlockTotal, filter);
+        rlBlockTotal = removeDataFromArrays(rlBlockTotal, filter);
+        rnBlockTotal = removeDataFromArrays(rnBlockTotal, filter);
+    }

    Matrix valid;
    if(transParams::applyCalib)
@@ -453,14 +460,14 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
    }
    dataInfno.findDefect = Matrix::New(findDefectData, 1, findDefectDataIndex);
    
-    // if(transParams::saveDebugInfomation)
-    // {
-    //     dataInfno.sn = snBlockTotal;
-    //     dataInfno.sl = slBlockTotal;
-    //     dataInfno.rn = rnBlockTotal;
-    //     dataInfno.rl = rlBlockTotal;
-    //     dataInfno.mp = mpBlockTotal;
-    // }
+    if(transParams::saveDebugInfomation)
+    {
+        dataInfno.sn = snBlockTotal;
+        dataInfno.sl = slBlockTotal;
+        dataInfno.rn = rnBlockTotal;
+        dataInfno.rl = rlBlockTotal;
+        dataInfno.mp = mpBlockTotal;
+    }

    tofDataTotal = removeDataFromArrays(tofDataTotal, valid);
    attDataTotal = removeDataFromArrays(attDataTotal, valid);
--- a/src/transmissionReconstruction/reconstruction/buildMatrix/FMM.cpp
+++ b/src/transmissionReconstruction/reconstruction/buildMatrix/FMM.cpp
@@ -13,7 +13,6 @@
 #include "Matrix.h"
 #include "Bresenham.h"
 #include "DGradient.h"
-#include "eikonal.h"


 using namespace Aurora;
--- a/src/transmissionReconstruction/reconstruction/reconstruction.cpp
+++ b/src/transmissionReconstruction/reconstruction/reconstruction.cpp
@@ -6,7 +6,6 @@
 #include "Matrix.h"
 #include "config/config.h"

-#include "CudaEnvInit.h"
 #include "log/notify.h"
 #include "transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h"
 #include "transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cuh"
@@ -136,10 +135,6 @@ namespace Recon {
        for (size_t i = 0; i < transParams::gpuSelectionList.getDataSize(); i++)
        {
            std::string msg;
-            if (!resetGPUDevice((int)transParams::gpuSelectionList[i],msg))
-            {
-                std::cerr<<msg<<std::endl;
-            }
        }
        std::vector<int> potentialMapDataDims = {1,1,1};
        for(size_t i=0; i<dims.getDataSize(); ++i)
@@ -195,9 +190,8 @@ namespace Recon {
        BuildMatrixResult buildMatrixR;
        for(int iter=1; iter<=numIter; ++iter)
        {
-            auto resDevice = res.toDeviceMatrix();
            //1200ms
-            buildMatrixR = buildMatrix(senderList.toDeviceMatrix(), receiverList.toDeviceMatrix(), resDevice, dims.toDeviceMatrix(), bentRecon && (iter!=1), potentialMap.toDeviceMatrix());
+            buildMatrixR = buildMatrix(senderList, receiverList, res, dims, bentRecon && (iter!=1), potentialMap);
            if(!data.isNull() && bentRecon && iter != numIter)
            {
                //与默认配置bentRecon不符，暂不实现todo
--- a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
+++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp
@@ -2,7 +2,6 @@

 #include "Function2D.h"
 #include "Matrix.h"
-#include "tval3gpu3d.h"

 #include "mkl_spblas.h"
 #include "mkl_solvers_ee.h"
@@ -21,99 +20,4 @@
 #include <sys/time.h>
 namespace Recon
 {
-    bool isEigsFinished = false;
-    bool isEigsStatus = false;
-
-    float eigs(Aurora::Sparse& aM)
-    {
-        size_t size = aM.getM();
-        if(size < aM.getN())
-        {
-            size = aM.getN();
-        }
-        Eigen::SparseMatrix<double> M(size,size);
-        Aurora::Matrix rows = aM.getRowVector();
-        Aurora::Matrix columns = aM.getColVector();
-        Aurora::Matrix values = aM.getValVector();
-        std::vector<Eigen::Triplet<double>> triplets(rows.getDataSize());
-        #pragma omp parallel for
-        for (int i = 0; i < rows.getDataSize(); ++i)
-        {
-            triplets[i] = Eigen::Triplet<double>(rows[i], columns[i], values[i]);
-        }
-        M.setFromTriplets(triplets.begin(), triplets.end());
-        float result;
-        Spectra::SparseGenMatProd<double> op(M);
-        Spectra::GenEigsSolver<Spectra::SparseGenMatProd<double>> eigs(op, 1, 6);
-        eigs.init();
-        int nconv = eigs.compute(Spectra::SortRule::LargestMagn);
-        Eigen::VectorXcd evalues;
-        if(eigs.info() == Spectra::CompInfo::Successful)
-        {
-            evalues = eigs.eigenvalues();
-            std::complex<double> complex = evalues[0];
-            result = complex.real();
-        }
-        isEigsFinished = true;
-        if (result> 1 + 1e-10)
-        {
-            isEigsStatus = true;
-        }
-        return result;
-    }
-
-    void checkAndScale(Aurora::Sparse& M, Aurora::Matrix& b,size_t n)
-    {
-        float s2 = 0.0;
-        if(!isEigsFinished)
-        {
-            s2 = eigs(M);
-            return;;
-        }
-
-        if (isEigsStatus)
-        {
-            b = b / sqrt(s2);
-            M.getValVector() = M.getValVector() / sqrt( s2);
-        }
-    }
-
-    Aurora::Matrix callTval3(Aurora::Sparse& M, Aurora::Matrix& b,const Aurora::Matrix& dims,int device, TVALOptions& opt)
-    {
-        checkAndScale(M,b,(size_t)Aurora::prod(dims).getScalar());     
-        int * yIdxs = new int[M.getColVector().getDataSize()];
-        std::copy(M.getColVector().getData(),M.getColVector().getData()+M.getColVector().getDataSize(),yIdxs);
-        int * xIdxs = new int[M.getRowVector().getDataSize()];
-        std::copy(M.getRowVector().getData(),M.getRowVector().getData()+M.getRowVector().getDataSize(),xIdxs);
-        Aurora::Matrix values = M.getValVector();
-        size_t cols = M.getM(), rows = M.getN();        
-        int nz = std::max(M.getColVector().getDataSize(),std::max(M.getRowVector().getDataSize(),M.getValVector().getDataSize()));        
-        sparse_matrix_t A;
-        sparse_matrix_t csrA;
-        mkl_sparse_s_create_coo(&A, sparse_index_base_t::SPARSE_INDEX_BASE_ZERO,  rows, cols,  nz,  yIdxs, xIdxs,values.getData());        
-        mkl_sparse_convert_csr(A, sparse_operation_t::SPARSE_OPERATION_NON_TRANSPOSE, &csrA);      
-        int n_rows,n_cols;
-        int *rows_start,*rows_end,*col_indx;
-        float * csrValues;
-        sparse_index_base_t index;
-        mkl_sparse_s_export_csr(csrA, &index, &n_rows, &n_cols, &rows_start, &rows_end, &col_indx, &csrValues);      
-        mkl_sparse_destroy(A);
-        delete [] xIdxs;
-        delete [] yIdxs;
-        
-        int *row_idx = new int[n_rows+1];
-        std::copy(rows_start,rows_start+n_rows,row_idx);
-        row_idx[n_rows] = rows_end[n_rows-1];
-        float* bData = b.getData();
-        // std::copy(b.getData(),b.getData()+b.getDataSize(),bData);
-        size_t bDims[3]={(size_t)b.getDimSize(0),(size_t)b.getDimSize(1),(size_t)b.getDimSize(2)};
-        size_t rdims[3] = {(size_t)dims[0], (size_t)dims[1], (size_t)dims[2]};
-        bool pagelocked = false;
-        auto result = TVALGPU( col_indx, row_idx,  csrValues, M.getM(), M.getN(), nz, bData, bDims, rdims, opt, device, false);
-        mkl_sparse_destroy(csrA);
-        delete [] row_idx;
-        //delete [] bData;
-        return Aurora::Matrix::fromRawData(result.data, result.dims[0],result.dims[1],result.dims[2]);
-        // return Aurora::Matrix();
-    }
 }
--- a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h
+++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h
@@ -1,17 +1,8 @@
 #ifndef __TVAL_H__
 #define __TVAL_H__
 #include "Sparse.h"
-#include "tvalstruct.h"

 namespace Recon {
-
-void checkAndScale(Aurora::Sparse& M, Aurora::Matrix& b,size_t n);
-
-float eigs(Aurora::Sparse& aM);
-
-Aurora::Matrix callTval3(Aurora::Sparse &M, Aurora::Matrix &b,
-                         const Aurora::Matrix &dims, int device,
-                         struct TVALOptions &options);
 }

 #endif // __TVAL_H__
--- a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.cpp
+++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.cpp
@@ -8,7 +8,6 @@
 #include "config/config.h"

 #include "transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h"
-#include "tvalstruct.h"

 namespace Recon
 {
@@ -26,25 +25,7 @@ namespace Recon
        if (Recon::transParams::name.empty()){
            Recon::transParams::name = "TVAL3";
        }
-        if (Recon::transParams::name == "TVAL3")
        {
-            //callTval3
-            TVALOptions opt;
-            opt.nonneg = solverOptions.nonNeg;
-            opt.bent = false;
-            opt.tol = 1E-10;
-            opt.maxit = niter;
-            opt.TVnorm = 2;
-            opt.disp = false;
-            opt.mu0 = solverOptions.TVAL3MU0;
-            opt.mu = solverOptions.TVAL3MU;
-            opt.beta = solverOptions.TVAL3Beta;
-            opt.beta0 = solverOptions.TVAL3Beta0;
-            int device = aDevice;
-            return callTval3(M, b, dims, device, opt);
-        }
-        //SART
-        else{
            //TODO:待实现，先实现默认的TVAL3
            return Aurora::Matrix();
        }
--- a/src/transmissionReconstruction/startTransmissionReconstruction.cpp
+++ b/src/transmissionReconstruction/startTransmissionReconstruction.cpp
@@ -25,9 +25,10 @@ TransmissionReconstructionResult Recon::startTransmissionReconstruction(const Au
                                     GeometryInfo& aGeomRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
                                     const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
 {
-    RECON_INFO("Start getTransmissionData.");
+    RECON_INFO("getTransmissionData start");
    auto transmissionData = getTransmissionData(aMotorPos, aMotoPosRef, aSlList, aSnList, aRlList, aRnList, aTemp, aTempRef, 
-                                                aGeom, aGeomRef, aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef);    
+                                                aGeom, aGeomRef, aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef);
+    RECON_INFO("getTransmissionData end");                                             
    Matrix dists = Recon::distanceBetweenTwoPoints(transmissionData.senderList, transmissionData.receiverList);
    Matrix sosRef = Recon::temperatureToSoundSpeed(transmissionData.waterTempList, "marczak");   
    //Recon::notifyProgress(17);
Author	SHA1	Message	Date
sunwen	f4c16d121b	feat: Add logic block log to performSignalProcessing in transmission reconstruction	2024-11-06 13:57:11 +08:00
kradchen	b859728562	feat: Add logic block log to performSignalProcessing in reflect reconstruction	2024-11-06 11:36:52 +08:00
kradchen	4bf226e4b2	feat: A CPU Performance Test tool branch, which remove cuda & outpout logic.	2024-11-04 10:55:21 +08:00
kradchen	f8d797da04	feat: Remove Cuda & other useless libs	2024-10-22 10:25:48 +08:00