12-22 commit for backup

2023-12-22 11:06:13 +08:00
parent 8ae5950e9f
commit 3c8268263c
13 changed files with 734 additions and 51 deletions
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -9,7 +9,10 @@ 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(BUILD_SHARED_LIBS ON)
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native")
+
+add_definitions(-DUSE_CUDA)
+# set(BUILD_SHARED_LIBS ON)
 if (${BUILD_SHARED_LIBS})
 set(MKL_INTERFACE_FULL intel_lp64)
 set(MKL_LINK static)
@@ -23,7 +26,7 @@ set(DCMTK_WITH_JPEG OFF)
 file(GLOB_RECURSE  cpp_files ./src/*.cpp)
 file(GLOB_RECURSE  cu_files ./src/*.cu)

-file(GLOB_RECURSE  header_files ./src/*.h)
+file(GLOB_RECURSE  header_files ./src/*.h ./src/*.cuh)

 if (${BUILD_SHARED_LIBS})
 set(cxx_files ./src/UR.cxx)
@@ -35,11 +38,10 @@ target_compile_options(UR PUBLIC $<TARGET_PROPERTY:MKL::MKL,INTERFACE_COMPILE_OP
 target_include_directories(UR PUBLIC $<TARGET_PROPERTY:MKL::MKL,INTERFACE_INCLUDE_DIRECTORIES>)
 target_link_libraries(UR PUBLIC $<LINK_ONLY:MKL::MKL>)
 target_link_libraries(UR PUBLIC OpenMP::OpenMP_CXX)
-target_compile_options(UR PUBLIC "-march=native")
 else()
 set(cxx_files ./src/main.cxx)
-add_executable(UR  ${cpp_files} ${cxx_files} ${header_files} ${Aurora_Source})
-target_compile_options(UR PUBLIC ${Aurora_Complie_Options} "-march=native")
+add_executable(UR  ${cpp_files} ${cu_files} ${cxx_files} ${header_files} ${Aurora_Source})
+target_compile_options(UR PUBLIC ${Aurora_Complie_Options} )
 endif()

 target_include_directories(UR PUBLIC ./src/)
@@ -71,7 +73,7 @@ set_target_properties(UR PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
 #                        -arch sm_75
 #                        >)
 target_compile_options(UR PRIVATE $<$<COMPILE_LANGUAGE:CUDA>:
-    -arch=sm_75
+    -arch=sm_75 --expt-extended-lambda 
 >)
 target_link_libraries(UR PRIVATE ${CUDA_RUNTIME_LIBRARY} CUDA::cufft CUDA::cudart)

@@ -82,23 +84,33 @@ 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 ./test/ ./src/ /usr/local/cuda/include)
+set_target_properties(UR_Test PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
+# target_compile_options(UR_Test PUBLIC ${Aurora_Complie_Options} "-march=native")
+target_include_directories(UR_Test PUBLIC ./src/)
 target_include_directories(UR_Test PUBLIC ${Aurora_INCLUDE_DIRS})
+target_include_directories(UR_Test PUBLIC ${DCMTK_INCLUDE_DIRS})
+target_include_directories(UR_Test PUBLIC ${Req_INCLUDES_DIRS})
 target_include_directories(UR_Test PUBLIC ${Parser_INCLUDE_DIRS})
 target_include_directories(UR_Test PUBLIC ${URDepends_INCLUDES_DIRS})
+target_include_directories(UR_Test PUBLIC ${DCMTK_INCLUDE_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 Req)
+target_link_libraries(UR_Test PUBLIC Parser)
 # 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)
+target_link_libraries(UR_Test PUBLIC dcmdata)

+target_compile_options(UR_Test PRIVATE $<$<COMPILE_LANGUAGE:CUDA>:
+    -arch=sm_75 --expt-extended-lambda 
+>)
+target_link_libraries(UR_Test PRIVATE ${CUDA_RUNTIME_LIBRARY} CUDA::cufft CUDA::cudart)
 gtest_discover_tests(UR_Test)
--- a/src/common/common.cpp
+++ b/src/common/common.cpp
@@ -63,6 +63,8 @@ namespace Recon
        auto coord_col1= aMCoordinates(Aurora::$,0).toMatrix();
        auto coord_col2Sub1= aMCoordinates(Aurora::$,1).toMatrix()-1;
        auto matrixSize_1 = aVMatrixSize.getData()[0];
+        printf("aMCoordinates size1:%zu\r\n",columns);
+
        switch (columns) {
            case 3:
            {
--- a/src/common/fileHelper.h
+++ b/src/common/fileHelper.h
@@ -6,8 +6,8 @@
 namespace Recon
 {
    const std::string DEFAULT_CONFIG_PATH = "/home/UR/ConfigFiles/";
-    const std::string DEFAULT_OUTPUT_PATH = "/home/UR/ReconResult/USCT_Result.mat";
-    const std::string DEFAULT_OUTPUT_FILENAME = "USCT_Result.mat";
+    const std::string DEFAULT_OUTPUT_PATH = "/home/UR/ReconResult/";
+    const std::string DEFAULT_OUTPUT_FILENAME = "USCT_Result1219.mat";

    std::string getPath(const std::string &aFullPath);
    bool endsWithMat(const std::string &aStr);
--- a/src/main.cxx
+++ b/src/main.cxx
@@ -3,11 +3,17 @@
 #include "config/config.h"
 #include "log/notify.h"
 #include "log/notify.h"
+#include <cstdio>
 #include <vector>
+#include "MatlabReader.h"
 #define EIGEN_USE_MKL_ALL

 #include "src/common/fileHelper.h"
 #include "startReconstructions.h"
+
+#include "transmissionReconstruction/detection/detection.h"
+#include "transmissionReconstruction/detection/detection.cuh"
+#include "startReconstructions.h"
 #include "log/log.h"
 /* 0 is data path.
    1 is dataRef path.
@@ -20,8 +26,8 @@ int main(int argc, char *argv[])
    int argNum = 5; 
    std::vector<std::string> args(argNum);
    args[0] = "";
-    args[1] = "";
-    args[2] = "";
+    args[1] = "/home/sun/20230418T145123/";
+    args[2] = "/home/sun/20230418T141000/";
    args[3] = Recon::DEFAULT_OUTPUT_PATH;
    args[4] = Recon::DEFAULT_CONFIG_PATH;
    argc = argc <= argNum? argc-1 : argNum;
@@ -29,30 +35,29 @@ int main(int argc, char *argv[])
    {
        args[i] = argv[i+1];
    }
-    if(args[0].empty())
-    {
-        printf("No recon id.");
-        return -1;
-    }
+    std::string outPutPath = args[3];
+    std::string directoryPath = outPutPath;
+    auto defaultLogger = getLogger("Main",outPutPath.data());
+    spdlog::set_default_logger(defaultLogger);
    std::string ReconID = args[0];
+    ReconID = ReconID=="none"?"":ReconID;
+    RECON_INFO("Read UR Args =====================");
+    RECON_INFO("ReconID:{0}",ReconID);
    if(args[1].empty())
    {
-        printf("No reconstruction data.");
+        RECON_INFO("No reconstruction data.");
        return -2;
    }
    std::string configPath = Recon::fixPathSlash(args[4]);
    Recon::initalizeConfig(configPath);
    if( args[2].empty() && Recon::transParams::runTransmissionReco)
    {
-        printf("Running transmission reconstruction, but no refrence data.");
+        RECON_INFO("Running transmission reconstruction, but no refrence data.");
        return -3;
    }
+RECON_INFO("configPath:{0}",configPath);

-    std::string outPutPath = args[3];
-    std::string directoryPath = outPutPath;
-    auto defaultLogger = getLogger("Main",outPutPath.data());
-    spdlog::set_default_logger(defaultLogger);
-
+    
    outPutPath = Recon::fixPathSlash(outPutPath);

    if(!Recon::isDirectory(directoryPath))
@@ -60,15 +65,18 @@ int main(int argc, char *argv[])
        RECON_INFO("Output directory is not valid.");
        return -4;
    }
-    
+    RECON_INFO("outPutPath:{0}",directoryPath);
    std::string dataPath = Recon::fixPathSlash(args[1]);
+RECON_INFO("dataPath:{0}",dataPath);
    std::string dataRefPath = Recon::fixPathSlash(args[2]);
-    RECON_INFO("start");
+    RECON_INFO("dataRefPath:{0}",dataRefPath);
+    RECON_INFO("UR Args End=======================");
+    RECON_INFO("UR Start");
    Recon::notifyStart(ReconID);
    int exitcode = Recon::startReconstructions(dataPath, dataRefPath, outPutPath);
    if (exitcode == 0)
    {
-        RECON_INFO("finish");
+        RECON_INFO("UR Finish");
        if (!Recon::notifyFinish()) {
            RECON_ERROR("Notify Finish failed!");
            return -100;
--- a/src/startReconstructions.cpp
+++ b/src/startReconstructions.cpp
@@ -33,7 +33,7 @@ using namespace Aurora;

 int Recon::startReconstructions(const std::string& aDataPath, const std::string& aDataRefPath, const std::string& aOutputPath)
 {
-    MatlabWriter writer(aOutputPath);
+    MatlabWriter writer(aOutputPath+"/result1220.mat");
    Parser dataParser(aDataPath);
    Parser refParser(aDataRefPath);
    Recon::DICOMExporter exporter(dataParser.getPatientData());
--- a/src/transmissionReconstruction/detection/detection.cpp
+++ b/src/transmissionReconstruction/detection/detection.cpp
@@ -2,6 +2,10 @@

 #include <algorithm>
 #include <cstddef>
+#include <cstdio>
+#include <cstdlib>
+#include <iostream>
+#include <ostream>
 #include <sys/types.h>

 #include "Function.h"
@@ -13,6 +17,8 @@
 #include "common/getMeasurementMetaData.h"
 #include "config/config.h"
 #include "calculateBankDetectAndHilbertTransformation.hpp"
+#include "log/log.h"
+#include "transmissionReconstruction/detection/detection.cuh"

 using namespace Aurora;
 namespace Recon {
@@ -84,7 +90,7 @@ namespace Recon {
        float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow)
    {
        auto sosOffset = calculateSOSOffset(sosBlock, expectedSOSWater, distBlock,  sampleRate);
-
+        
        auto calcResult = calculateStarEndSearchPosition(distBlock, minSpeedOfSound, maxSpeedOfSound, sampleRate,AscanBlock.getDimSize(0), sosOffset, startOffset, segmentLenOffset);

        auto AscanBlockProcessed = zeros(AscanBlock.getDimSize(0),AscanBlock.getDimSize(1));
@@ -302,7 +308,7 @@ namespace Recon {
            max(c($, i).toMatrix(), All, rowID, colID);
            shiftInSamples[i] = -maxlag + colID + rowID;
        }
-        if (useTimeWindowing) {
+                if (useTimeWindowing) {
            shiftInSamples = shiftInSamples - diffStartSearch;
        }
        auto tof = shiftInSamples / sampleRate + distBlock / sosWaterBlock;
@@ -339,10 +345,10 @@ namespace Recon {
        TimeWindowResult timeResult2;
        timeResult2.AscanBlockProcessed = AscanRefBlock;
        if (useTimeWindowing == 1) {
-            timeResult1 = applyTimeWindowing(
-                AscanBlock, sampleRate, distBlock, sosWaterBlock,
-                expectedSOSWater, offsetElectronicSamples, detectionWindowSOS,
-                minSpeedOfSound, maxSpeedOfSound, gaussWindow);
+                        timeResult1 = applyTimeWindowing(
+                    AscanBlock, sampleRate, distBlock, sosWaterBlock,
+                    expectedSOSWater, offsetElectronicSamples, detectionWindowSOS,
+                    minSpeedOfSound, maxSpeedOfSound, gaussWindow);
            timeResult2 = applyTimeWindowing(
                AscanRefBlock, sampleRate, distRefBlock, sosWaterRefBlock,
                expectedSOSWater, offsetElectronicSamples, detectionWindowSOS,
@@ -424,24 +430,29 @@ namespace Recon {
        auto _sosWaterBlock = temperatureToSoundSpeed(sosWaterBlock, "marczak");
        auto _sosWaterRefBlock = temperatureToSoundSpeed(sosWaterRefBlock, "marczak");
        switch (Recon::transParams::version) {
-        case 1: {
-            return detectTofAndAttMex(
-                AscanBlock, AscanRefBlock, distBlock, distRefBlock,
-                _sosWaterBlock, _sosWaterRefBlock, Recon::transParams::resampleFactor, Recon::transParams::nThreads,
-                expectedSOSWater, Recon::transParams::useTimeWindowing,
-                Recon::transParams::aScanReconstructionFrequency, Recon::transParams::detectionWindowATT,
-                Recon::transParams::offsetElectronic, Recon::transParams::detectionWindowSOS, Recon::transParams::minSpeedOfSound,
-                Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow);
-        }
-        case 2:
+        // case 1: {
+        //     return detectTofAndAttMex(
+        //         AscanBlock, AscanRefBlock, distBlock, distRefBlock,
+        //         _sosWaterBlock, _sosWaterRefBlock, Recon::transParams::resampleFactor, Recon::transParams::nThreads,
+        //         expectedSOSWater, Recon::transParams::useTimeWindowing,
+        //         Recon::transParams::aScanReconstructionFrequency, Recon::transParams::detectionWindowATT,
+        //         Recon::transParams::offsetElectronic, Recon::transParams::detectionWindowSOS, Recon::transParams::minSpeedOfSound,
+        //         Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow);
+        // }
+        // case 2:
        default:
-            return detectTofAndAtt(
-                AscanBlock, AscanRefBlock, distBlock, distRefBlock,
-                _sosWaterBlock, _sosWaterRefBlock, Recon::transParams::resampleFactor, Recon::transParams::nThreads,
+            auto r =  detectTofAndAtt(
+                AscanBlock.toDeviceMatrix(), AscanRefBlock.toDeviceMatrix(), distBlock.toDeviceMatrix(), distRefBlock.toDeviceMatrix(),
+                _sosWaterBlock.toDeviceMatrix(), _sosWaterRefBlock.toDeviceMatrix(), Recon::transParams::resampleFactor, Recon::transParams::nThreads,
                expectedSOSWater, Recon::transParams::useTimeWindowing,
                Recon::transParams::aScanReconstructionFrequency, Recon::transParams::detectionWindowATT,
                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();
+            return ret;
        }
    }
 }
--- a/src/transmissionReconstruction/detection/detection.cu
+++ b/src/transmissionReconstruction/detection/detection.cu
@@ -0,0 +1,323 @@
+#include "CudaMatrix.h"
+#include "detection.cuh"
+#include <math.h>
+#include <thrust/transform.h>
+#include <thrust/execution_policy.h>
+#include "Function1D.cuh"
+#include "Function2D.cuh"
+#include "Function3D.h"
+
+__global__ void copyMatrixKernel(float* aStartPos1, float* aEndPos1,
+                                float* aStartPos2, float* aEndPos2,
+                                float* aIn1, float* aOut1, 
+                                float* aIn2, float* aOut2, 
+                                unsigned int aColElements){
+    unsigned int i = blockIdx.x;
+    unsigned int base_offset = i * aColElements;
+    unsigned int j = aStartPos1[i]-1+threadIdx.x;
+    for (int offset = 0; j+offset<aEndPos1[i]; offset+=blockDim.x) {
+        aOut1[j + base_offset + offset] = aIn1[j + base_offset + offset];
+    }
+    j = aStartPos2[i]-1+threadIdx.x;
+    for (int offset = 0; j+offset<aEndPos2[i]; offset+=blockDim.x) {
+        aOut2[j + base_offset + offset] = aIn2[j + base_offset + offset];
+    }
+
+}
+
+__global__ void setPosKernel(float* aStartPos1, float* aEndPos1,
+                                float* aStartPos2, float* aEndPos2,
+                                float* tof, float* tof2, 
+                                float* sizeAscan, 
+                                int sampleRate,
+                                float offsetElectronicSamples,
+                                int detectionWindowATT)
+{
+    unsigned int i = blockIdx.x;
+    aStartPos1[i] = floor(fmaxf(tof[i] * sampleRate + offsetElectronicSamples, (float)1.0));
+    aEndPos1[i] = ceilf(fminf(sizeAscan[0], tof[i] * sampleRate + offsetElectronicSamples +
+                              detectionWindowATT));
+    aStartPos2[i] = floorf(fmaxf(tof2[i], 1.0f));
+    aEndPos2[i] = ceilf(fminf(sizeAscan[1], tof2[i] + detectionWindowATT));
+}
+
+
+CudaMatrix Recon::calculateAttenuationCuda(const CudaMatrix &ascans,
+                                        const CudaMatrix &startPos,
+                                        const CudaMatrix &endPos,
+                                        const CudaMatrix &ascansRef,
+                                        const CudaMatrix &startPosRef,
+                                        const CudaMatrix &endPosRef){
+    auto ascans2 = Aurora::zerosCuda(ascans.getDimSize(0), ascans.getDimSize(1));
+    auto ascansRef2 = Aurora::zerosCuda(ascansRef.getDimSize(0), ascansRef.getDimSize(1));
+    copyMatrixKernel<<<ascans2.getDimSize(1),256>>>(startPos.getData(),endPos.getData(),
+                                                    startPosRef.getData(), endPosRef.getData(),
+                                                    ascans.getData(), ascans2.getData(),
+                                                    ascansRef.getData(), ascansRef2.getData(),
+                                                    ascans2.getDimSize(0));
+
+    auto pulseEnergy = Aurora::sum(ascans2^2);
+    auto pulseEnergyEmpty = Aurora::sum(ascansRef2^2);
+
+    return Aurora::log(pulseEnergyEmpty/pulseEnergy);
+}
+
+CudaMatrix
+Recon::detectAttVectorizedCuda(const CudaMatrix &Ascan, const CudaMatrix &AscanRef,
+                    const CudaMatrix &distRef,
+                    const CudaMatrix &sosWaterRef,
+                    const CudaMatrix &tof, int aScanReconstructionFrequency,
+                    float offsetElectronic, int detectionWindowATT) {
+    auto sizeAscan = size(Ascan);
+
+    auto sampleRate = aScanReconstructionFrequency;
+    float offsetElectronicSamples = offsetElectronic * sampleRate;
+
+    auto envelopeOfAScan = Aurora::abs(Aurora::hilbert(Ascan));
+    auto envelopeOfReferenceAScan = Aurora::abs(hilbert(AscanRef));
+    auto tof2 = (distRef / sosWaterRef) * sampleRate + offsetElectronicSamples;
+
+    auto startPos = zerosCuda(Ascan.getDimSize(1), 1);
+    auto endPos = zerosCuda(Ascan.getDimSize(1), 1);
+    auto startPosRef = zerosCuda(Ascan.getDimSize(1), 1);
+    auto endPosRef = zerosCuda(Ascan.getDimSize(1), 1);
+
+    setPosKernel<<<Ascan.getDimSize(1),1>>>(startPos.getData(),
+                                endPos.getData(), startPosRef.getData(), endPosRef.getData(),
+                                tof.getData(), tof2.getData(), sizeAscan.getData(), sampleRate,
+                                offsetElectronicSamples,detectionWindowATT);
+    return calculateAttenuationCuda(envelopeOfAScan, startPos, endPos,
+                                envelopeOfReferenceAScan, startPosRef,
+                                endPosRef);
+}
+
+CudaMatrix Recon::calculateSOSOffset(const CudaMatrix &sosBlock, float referenceSOS, const CudaMatrix &distBlock, float sampleRate){
+    auto offset = (distBlock / sosBlock - distBlock / referenceSOS) * sampleRate;
+    return offset;
+}
+
+Recon::SearchPositionC Recon::calculateStarEndSearchPosition(const CudaMatrix &aVDistBlock,
+                                float minSpeedOfSound, float maxSpeedOfSound,
+                                float sampleRate, float maxSample,
+                                const CudaMatrix &aVSosOffsetBlock,
+                                float startOffset, float segmentLenOffset)
+{
+        auto startSearch = floor((aVDistBlock / maxSpeedOfSound)*sampleRate+aVSosOffsetBlock+startOffset);
+        auto lambda = [=] __device__ (const float& v){
+            float a = (uint)(v);
+            if(a<1) return 1.0f;
+            else if(a>maxSample) return maxSample;
+            return a;
+        };
+        thrust::transform(thrust::device, startSearch.getData(),startSearch.getData()+startSearch.getDataSize(),
+                        startSearch.getData(),lambda);
+
+        auto endSearch = ceil(aVDistBlock/minSpeedOfSound*sampleRate+aVSosOffsetBlock+startOffset+segmentLenOffset);
+
+        thrust::transform(thrust::device, endSearch.getData(),endSearch.getData()+endSearch.getDataSize(),
+                        endSearch.getData(),lambda);
+        Recon::SearchPositionC result;
+        result.startSearch = startSearch;
+        result.endSearch = endSearch;
+        return result;
+}
+
+__global__ void copyMatrixKernel(float* aStartPos1, float* aEndPos1,
+                                float* aIn1, float* aOut1, 
+                                unsigned int aColElements){
+    unsigned int i = blockIdx.x;
+    unsigned int base_offset = i * aColElements;
+    unsigned int j = aStartPos1[i]-1+threadIdx.x;
+    for (int offset = 0; j+offset<aEndPos1[i]; offset+=blockDim.x) {
+        aOut1[j + base_offset + offset] = aIn1[j + base_offset + offset];
+    }
+}
+
+__device__ float getGuassValue(int aWidth,int aIndex){
+    float v = (-5.0f+(10.0f/(float)(aWidth-1))*aIndex);
+    return exp(-0.1f*(v*v));
+}
+
+__global__ void guassWindowKernel(float* aStartSearch,float* aEndSearch,
+                                float* aIn1, float* aOut1, float* aExpectedPosWater,
+                                unsigned int aColElements){
+    __shared__ unsigned int windowWidth;
+    __shared__ unsigned int base_offset;
+    __shared__ unsigned int j;
+    __shared__ unsigned int begin;
+    __shared__ unsigned int end ;
+
+    if(threadIdx.x == 0){
+        windowWidth = aEndSearch[blockIdx.x]-aStartSearch[blockIdx.x];
+        base_offset = blockIdx.x * aColElements;
+        j = aStartSearch[blockIdx.x]-1+threadIdx.x;
+        begin = roundf(aExpectedPosWater[blockIdx.x])- roundf(windowWidth/2)-1;
+        end = roundf(aExpectedPosWater[blockIdx.x])-roundf(windowWidth/2)+windowWidth-1;
+    }
+
+    for (int offset = 0; j+offset<aEndSearch[blockIdx.x]; offset+=blockDim.x) {
+        unsigned int threadAccIdx = j + offset;
+        if(threadAccIdx>=begin && threadAccIdx<end){
+            aOut1[threadAccIdx + base_offset] = aIn1[threadAccIdx + base_offset]*getGuassValue(windowWidth, threadAccIdx - begin);
+        }
+        else{
+            aOut1[threadAccIdx + base_offset] = aIn1[threadAccIdx + base_offset];
+        }
+    }
+}
+
+Recon::TimeWindowResultC Recon::applyTimeWindowing(const Aurora::CudaMatrix &AscanBlock, float sampleRate,
+        const Aurora::CudaMatrix &distBlock, const Aurora::CudaMatrix &sosBlock,
+        float expectedSOSWater, float startOffset, float segmentLenOffset,
+        float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow)
+    {
+        auto sosOffset = calculateSOSOffset(sosBlock, expectedSOSWater, distBlock,  sampleRate);
+
+        auto calcResult = calculateStarEndSearchPosition(distBlock, minSpeedOfSound, maxSpeedOfSound, sampleRate,AscanBlock.getDimSize(0), sosOffset, startOffset, segmentLenOffset);
+
+        auto AscanBlockProcessed = zerosCuda(AscanBlock.getDimSize(0),AscanBlock.getDimSize(1));
+
+        if(gaussWindow)
+        {
+            auto expectedPosWater = (distBlock / expectedSOSWater) * sampleRate + startOffset;
+            guassWindowKernel<<<AscanBlock.getDimSize(1),256>>>(calcResult.startSearch.getData(),
+                            calcResult.endSearch.getData(), AscanBlock.getData(), AscanBlockProcessed.getData(),
+                            expectedPosWater.getData(), AscanBlock.getDimSize(0));
+        }
+        else{
+            copyMatrixKernel<<<AscanBlock.getDimSize(1),256>>>(calcResult.startSearch.getData(),
+                            calcResult.endSearch.getData(),AscanBlock.getData(),AscanBlockProcessed.getData(),
+                            AscanBlock.getDimSize(0));
+        }
+        Recon::TimeWindowResultC result;
+        result.startSearch = calcResult.startSearch;
+        result.AscanBlockProcessed = AscanBlockProcessed;
+        return result;
+    }
+
+__device__ struct KVPair{
+    __device__ KVPair(){}
+    __device__ KVPair(int aIndex,float aValue):index(aIndex),value(aValue){};
+    __device__ KVPair(const KVPair& other){
+        index = other.index;
+        value = other.value;
+    }
+    int index;
+    float value;
+};
+
+__device__ KVPair maxKV(const KVPair& x, const KVPair& y){
+    return x.value>y.value ? x:y;
+};
+
+__global__ void findMaxIndexKernel(float* aInputData, unsigned int aColSize,float* aOutput, int aMaxlag){
+    //确定每个thread的index
+    unsigned int idx = blockIdx.x * aColSize + threadIdx.x;
+    __shared__ KVPair shared_data[256];
+    // 每个线程加载一个元素到共享内存
+    shared_data[threadIdx.x] = (threadIdx.x< aColSize) ? KVPair(threadIdx.x,aInputData[idx]) : KVPair(0,-FLT_MAX);
+    __syncthreads();
+    // 每个线程规约自己的分段，将每个blockDim.x的值规约到数组最前面一段
+    for (int offset = blockDim.x; offset<aColSize; offset+=blockDim.x) {
+        if(threadIdx.x + offset<aColSize){
+            shared_data[threadIdx.x] = maxKV(shared_data[threadIdx.x], KVPair(threadIdx.x+offset, aInputData[idx + offset]));
+        }
+        __syncthreads();
+    }
+    // 规约最前面一段
+    for (int i = blockDim.x/2; i >0; i>>=1) {
+
+        if (threadIdx.x < i) {
+            shared_data[threadIdx.x] = maxKV(shared_data[threadIdx.x], shared_data[i + threadIdx.x]);
+        }
+        __syncthreads();
+    }
+
+    // 第一个线程存储每个分段的最大值到全局内存
+    if (threadIdx.x == 0) {
+        aOutput[blockIdx.x] = shared_data[0].index - aMaxlag;
+    }
+}
+int nextpow2(unsigned int value){
+        unsigned int compareValue = 2;
+        int result = 1;    
+        while(compareValue<=value){
+            compareValue=compareValue<<1;
+            result++;
+        }
+        return result;
+}
+Recon::DetectResultC Recon::detectTofAndAtt(
+        const Aurora::CudaMatrix &AscanBlock, const Aurora::CudaMatrix &AscanRefBlock,
+        const Aurora::CudaMatrix &distBlock, const Aurora::CudaMatrix &distRefBlock,
+        const Aurora::CudaMatrix &sosWaterBlock,
+        const Aurora::CudaMatrix &sosWaterRefBlock,
+        int resampleFactor,int nthreads, float expectedSOSWater,
+        int useTimeWindowing, int aScanReconstructionFrequency,int detectionWindowATT,
+        float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound,
+        float maxSpeedOfSound, bool gaussWindow)
+    {
+        auto sampleRate = aScanReconstructionFrequency;
+        float offsetElectronicSamples = offsetElectronic * sampleRate;
+        CudaMatrix diffStartSearch;
+        Recon::TimeWindowResultC timeResult1;
+        timeResult1.AscanBlockProcessed = AscanBlock;
+        Recon::TimeWindowResultC timeResult2;
+        timeResult2.AscanBlockProcessed = AscanRefBlock;
+        if (useTimeWindowing == 1) {
+            timeResult1 = applyTimeWindowing(
+                AscanBlock, sampleRate, distBlock, sosWaterBlock,
+                expectedSOSWater, offsetElectronicSamples, detectionWindowSOS,
+                minSpeedOfSound, maxSpeedOfSound, gaussWindow);
+            timeResult2 = applyTimeWindowing(
+                AscanRefBlock, sampleRate, distRefBlock, sosWaterRefBlock,
+                expectedSOSWater, offsetElectronicSamples, detectionWindowSOS,
+                minSpeedOfSound, maxSpeedOfSound, gaussWindow);
+
+            diffStartSearch = timeResult1.startSearch - timeResult2.startSearch;
+        }
+        auto _AscanBlock = timeResult1.AscanBlockProcessed;
+        auto _AscanRefBlock = timeResult2.AscanBlockProcessed;
+        int maxlag = 0;
+        CudaMatrix c, c1;
+        {
+            auto m = std::max(size(_AscanBlock, 1), size(_AscanRefBlock, 1));
+            maxlag =
+                std::max(size(_AscanBlock, 1), size(_AscanRefBlock, 1)) - 1;
+            auto mxl = std::min(maxlag, m - 1);
+            auto ceilLog2 = nextpow2(2 * m - 1);
+            auto m2 = pow(2, ceilLog2);
+            {
+                CudaMatrix c_1_2;
+                {
+                    CudaMatrix c_1_1;
+                    {
+                        auto x = fft(_AscanBlock, m2);
+                        auto y = fft(_AscanRefBlock, m2);
+                        c_1_1 = x * conj(y);
+                    }
+                    c_1_2 = ifft(c_1_1);
+                } 
+                c1 = real(c_1_2);
+            }
+            c = zerosCuda(mxl + mxl + 1, c1.getDimSize(1));
+            c.setBlock(0, 0, mxl-1, c1.block(0, m2-mxl,m2-1));
+            c.setBlock(0, mxl, mxl*2, c1.block(0, 0, mxl));
+        }
+        auto shiftInSamples = zerosCuda(1, c1.getDimSize(1));
+        findMaxIndexKernel<<<c.getDimSize(1),256>>>(c.getData(),c.getDimSize(0),shiftInSamples.getData(),maxlag);
+        if (useTimeWindowing) {
+            shiftInSamples = shiftInSamples - diffStartSearch;
+        }
+        auto tof = shiftInSamples / sampleRate + distBlock / sosWaterBlock;
+        auto sosValue = distBlock / tof;
+        Recon::DetectResultC result;
+        result.sosValue = sosValue;
+        auto tofRel = tof - distBlock / sosWaterBlock;
+        result.att = detectAttVectorizedCuda(
+            _AscanBlock, _AscanRefBlock, distRefBlock, sosWaterRefBlock,
+            tof, aScanReconstructionFrequency, offsetElectronic,
+            detectionWindowATT);
+        result.tof = tofRel;
+        return result;
+    }
--- a/src/transmissionReconstruction/detection/detection.cuh
+++ b/src/transmissionReconstruction/detection/detection.cuh
@@ -0,0 +1,53 @@
+#ifndef __DETECTION_CUH__
+#define __DETECTION_CUH__
+#include "CudaMatrix.h"
+using namespace Aurora;
+namespace Recon {
+struct SearchPositionC {
+  Aurora::CudaMatrix startSearch;
+  Aurora::CudaMatrix endSearch;
+};
+struct DetectResultC {
+  Aurora::CudaMatrix tof;
+  Aurora::CudaMatrix sosValue;
+  Aurora::CudaMatrix att;
+};
+struct TimeWindowResultC {
+  Aurora::CudaMatrix startSearch;
+  Aurora::CudaMatrix AscanBlockProcessed;
+};
+    CudaMatrix calculateAttenuationCuda(const CudaMatrix &ascans,
+                                        const CudaMatrix &startPos,
+                                        const CudaMatrix &endPos,
+                                        const CudaMatrix &ascansRef,
+                                        const CudaMatrix &startPosRef,
+                                        const CudaMatrix &endPosRef);
+    CudaMatrix
+    detectAttVectorizedCuda(const CudaMatrix &Ascan, const CudaMatrix &AscanRef,
+                        const CudaMatrix &distRef,
+                        const CudaMatrix &sosWaterRef,
+                        const CudaMatrix &tof, int aScanReconstructionFrequency,
+                        float offsetElectronic, int detectionWindowATT);
+    
+    TimeWindowResultC applyTimeWindowing(const Aurora::CudaMatrix &AscanBlock, float sampleRate,
+        const Aurora::CudaMatrix &distBlock, const Aurora::CudaMatrix &sosBlock,
+        float expectedSOSWater, float startOffset, float segmentLenOffset,
+        float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow);
+    SearchPositionC calculateStarEndSearchPosition(const CudaMatrix &aVDistBlock,
+                                float minSpeedOfSound, float maxSpeedOfSound,
+                                float sampleRate, float maxSample,
+                                const CudaMatrix &aVSosOffsetBlock,
+                                float startOffset, float segmentLenOffset);
+    CudaMatrix calculateSOSOffset(const CudaMatrix &sosBlock, float referenceSOS, const CudaMatrix &distBlock, float sampleRate);
+    DetectResultC detectTofAndAtt(
+        const Aurora::CudaMatrix &AscanBlock, const Aurora::CudaMatrix &AscanRefBlock,
+        const Aurora::CudaMatrix &distBlock, const Aurora::CudaMatrix &distRefBlock,
+        const Aurora::CudaMatrix &sosWaterBlock,
+        const Aurora::CudaMatrix &sosWaterRefBlock,
+        int resampleFactor,int nthreads, float expectedSOSWater,
+        int useTimeWindowing, int aScanReconstructionFrequency,int detectionWindowATT,
+        float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound,
+        float maxSpeedOfSound, bool gaussWindow);
+};
+
+#endif // __DETECTION_H__
--- a/src/transmissionReconstruction/reconstruction/buildMatrix/Bresenham.cpp
+++ b/src/transmissionReconstruction/reconstruction/buildMatrix/Bresenham.cpp
@@ -1,4 +1,5 @@
 #include "Bresenham.h"
+#include <cstdio>
 #include <math.h> 
 #include <stdlib.h>
 #include <vector>
@@ -109,10 +110,85 @@ float* b3dMex(float *startPoint, float *endPoint,size_t &outputSize ) {
    outputSize = points.size();
    return output;
 }
+float * b3dMexDouble2(float *startPoint,float *endPoint, size_t& outputSize) {
+
+    std::vector<float> points; 
+    
+    float i, l, m, n, x_inc, y_inc, z_inc,err_0, err_1, err_2, dx2, dy2, dz2,d, x,y,z;
+           
+    x = startPoint[0];
+    y = startPoint[1];
+    z = startPoint[2];
+    
+    float dx = endPoint[0] - x;
+    float dy = endPoint[1] - y;
+    float dz = endPoint[2] - z;
+    
+    x_inc = (dx < 0) ? -1 : 1;
+    l = abs(dx);
+    y_inc = (dy < 0) ? -1 : 1;
+    m = abs(dy);
+    z_inc = (dz < 0) ? -1 : 1;
+    n = abs(dz);
+    dx2 = l*2;// << 1;
+    dy2 = m*2;// << 1;
+    dz2 = n*2;// << 1;
+    float maxV;
+    if ((l >= m) && (l >= n)){
+        maxV = l;
+        d = dx2;
+        dx2 = 0;
+    }
+    else if((m >= l) && (m >= n)){
+        maxV = m;
+        d = dy2;
+        dy2 = 0;
+    }
+    else{
+        maxV = n;
+        d = dz2;
+        dz2 = 0;
+    }
+    err_0 = dx2==0?0:(dx2 - maxV);
+    err_1 = dy2==0?0:(dy2 - maxV);
+    err_2 = dz2==0?0:(dz2 - maxV);
+
+    for (i = 0; i < maxV; i++) {
+
+        points.push_back(x);
+        points.push_back(y);
+        points.push_back(z);
+        if (err_0 > 0) {
+          x += x_inc;
+          err_0 -= d;
+        }
+        if (err_1 > 0) {
+          y += y_inc;
+          err_1 -= d;
+        }
+        if (err_2 > 0) {
+          z += z_inc;
+          err_2 -= d;
+        }
+        err_0 += dx2;
+        err_1 += dy2;
+        err_2 += dz2;
+        if (err_0 == 0)x += x_inc;
+        if (err_1 == 0)y += y_inc;
+        if (err_2 == 0)z += z_inc;
+
+    }
+    float* output = new float[points.size()];
+    std::copy(points.begin(), points.end(), output);
+    outputSize = points.size();
+    return output;
+}

 Aurora::Matrix b3dMexDouble(const Aurora::Matrix& StartPt, const Aurora::Matrix& endPts)
 {
    size_t size = 0;
+    size_t size1 = 0;
+
    auto temp  =
        Recon::b3dMexDouble(StartPt.getData(), endPts.getData(), size);
    return Aurora::Matrix::fromRawData(temp, 3, size / 3);
@@ -215,5 +291,7 @@ float * b3dMexDouble(float *startPoint,float *endPoint, size_t& outputSize) {
    outputSize = points.size();
    return output;
 }
+
+
 }

--- a/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp
+++ b/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp
@@ -36,6 +36,7 @@ namespace Recon
        auto path = transpose(b3dMexDouble(startPt, endPt));
        result.pathLen = path.getDimSize(0);
        //uint32(path)
+        //75,3
        result.path = round(path);
        result.weighting = getPixelLengthApproximation(startPt, endPt, res, result.pathLen);
        return result;
@@ -270,12 +271,14 @@ namespace Recon
        for (size_t ray = 0; ray < nTotalRays; ray++)
        {
            rayCount++;
+            //1, 3, 1
            auto startPt = transpose(senderList($,ray).toMatrix());
            auto endPt = transpose(receiverList($,ray).toMatrix());
            if (!BENT){
                if (Recon::transParams::bresenham)
                {
                    auto trace = traceStraightRayBresenham(startPt, endPt, res);
+                    // 75, 3
                    path = trace.path;
                    weighting = trace.weighting;
                    pathLenDisc = trace.pathLen;
@@ -295,6 +298,7 @@ namespace Recon
                    break;
                }
            }
+            printf("dims size1:%zu\r\n",dims.getDimSize(1));
            Matrix linearIndices =
                dims.getDimSize(1) == 2
                    ? convertToLinearIndices(dims, path.block(1, 0, 1))
@@ -302,9 +306,9 @@ namespace Recon
            linearIndices = linearIndices-1;
            linearIndices.forceReshape(1, linearIndices.getDataSize(), 1);
            if (Recon::transParams::saveDebugInfomation){
-                for (size_t i = 0; i < linearIndices.getDataSize(); i++)
+                for (size_t n = 0; n < linearIndices.getDataSize(); n++)
                {
-                    result.hitmap[linearIndices[i]]+=1;
+                    result.hitmap[linearIndices[n]]+=1;
                }
            }
            // printf("Progress: %f (%zu of %zu)\r\n",(float)rayCount*100/(float)nTotalRays,rayCount,nTotalRays);
--- a/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cu
+++ b/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cu
@@ -0,0 +1,181 @@
+#include "CudaMatrix.h"
+#include <thrust/reduce.h>
+#include <thrust/execution_policy.h>
+#include <cuda_runtime.h>
+
+#include <math.h>
+namespace  {
+    const int THREADS_PER_BLOCK = 256;
+}
+
+
+
+__device__ float getPixelLengthKernel(float* aStartPt, float*aEndPt, float* aRes, int aPathLenDisc)
+{
+    float3 temp;
+    temp.x = aEndPt[0]- aStartPt[0];
+    temp.y = aEndPt[1]- aStartPt[1];
+    temp.z = aEndPt[2]- aStartPt[2];
+    temp.x *= aRes[0];
+    temp.y *= aRes[1];
+    temp.z *= aRes[2];
+    temp.x *= temp.x;
+    temp.y *= temp.y;
+    temp.z *= temp.x;
+    return sqrtf(temp.x+temp.y+temp.z)/aPathLenDisc;
+}
+
+__global__ struct b3dStruct{
+    float max, x_inc, y_inc, z_inc,err_0, err_1, err_2, dx2, dy2, dz2, d, x, y, z, weight;
+};
+__device__ float * getRayBLen(float *startPoint, float *endPoint, struct b3dStruct* aOut) {
+       
+    aOut->x = startPoint[0];
+    aOut->y = startPoint[1];
+    aOut->z = startPoint[2];
+    
+    float dx = endPoint[0] - aOut->x;
+    float dy = endPoint[1] - aOut->y;
+    float dz = endPoint[2] - aOut->z;
+    
+    aOut->x_inc = (dx < 0) ? -1 : 1;
+    float l = abs(dx);
+    aOut->y_inc = (dy < 0) ? -1 : 1;
+    float m = abs(dy);
+    aOut->z_inc = (dz < 0) ? -1 : 1;
+    float n = abs(dz);
+    aOut->dx2 = l*2;// << 1;
+    aOut->dy2 = m*2;// << 1;
+    aOut->dz2 = n*2;// << 1;
+    if ((l >= m) && (l >= n)){
+        aOut->max = l;
+        aOut->d = aOut->dx2;
+        aOut->dx2 = 0;
+    }
+    else if((m >= l) && (m >= n)){
+        aOut->max = m;
+        aOut->d = aOut->dy2;
+        aOut->dy2 = 0;
+    }
+    else{
+        aOut->max = n;
+        aOut->d = aOut->dz2;
+        aOut->dz2 = 0;
+    }
+    aOut->err_0 = aOut->dx2==0?0:(aOut->dx2 - aOut->max);
+    aOut->err_1 = aOut->dy2==0?0:(aOut->dy2 - aOut->max);
+    aOut->err_2 = aOut->dz2==0?0:(aOut->dz2 - aOut->max);
+}
+
+/**
+ * @brief 计算traceStraightRayBresenham的pathlen和weight的核函数
+ * 
+ * @param aStratPt 
+ * @param aEndPt 
+ * @param aSize 
+ * @param aRes 
+ * @param aOuts 
+ * @return _
+ */
+__global__ void calcLenWeightKernel(float* aStratPt, float* aEndPt, int aSize,float* aRes, b3dStruct* aOuts ){
+    unsigned int idx = blockIdx.x*blockDim.x + threadIdx.x;
+    if (idx > aSize)return;
+    getRayBLen(aStratPt + idx*3, aEndPt + idx*3, aOuts+idx);
+    aOuts[idx].weight = getPixelLengthKernel(aStratPt + idx*3, aEndPt + idx*3, aRes + idx*3,aOuts[idx].max);
+};
+
+__global__ void calcPathOffsetKernel(b3dStruct* aStructs ,float* aPath, int aSize){
+    aPath[0] = aStructs[0].max;
+    for (int i=1; i<aSize; i++) {
+        aPath[i] =aStructs[i].max+aPath[i]-1; 
+    }
+}
+
+__global__ void convertToLinearIndices(const float3& aDims, const float3& aPath, float* aOut){
+    float matrixSize_1 = aDims.x;
+    float matrixSize_2 = aDims.y;
+    float coord_col1 = aPath.x;
+    float coord_col2Sub1 = aPath.y - 1;
+    float coord_col3Sub1 = aPath.z - 1;
+    *aOut = coord_col1 + coord_col2Sub1 * matrixSize_1 + coord_col3Sub1 * matrixSize_1 * matrixSize_2 -1;
+}
+
+__global__ void calcRayBLPath(struct b3dStruct* aInstruct,
+            float* dims, unsigned int * aOffset, float* aIOut,
+            float* aJOut, float* aSOut, unsigned int aSize) { 
+    unsigned int idx = blockIdx.x*blockDim.x + threadIdx.x;
+    if (idx > aSize)return;
+        b3dStruct* in = aInstruct+idx;
+float max, x_inc, y_inc, z_inc,err_0, err_1, err_2, dx2, dy2, dz2,d, x,y,z;
+    x = in->x;
+    y = in->y;
+    z = in->z;
+    x_inc = in->x_inc;
+    y_inc = in->y_inc;
+    z_inc = in->z_inc;
+    err_0 = in->err_0;
+    err_1 = in->err_1;
+    err_2 = in->err_2;
+    max = in->max;
+
+    float3 dims3{dims[0],dims[1],dims[2]};
+    float* output_ptr = aOutput + aOffset[idx];
+    for (float i = 0; i < max; i++) {
+        convertToLinearIndices(dims3, {roundf(x),roundf(y),roundf(z)}, output_ptr);
+        output_ptr ++;
+        if (err_0 > 0) {
+          x += x_inc;
+          err_0 -= d;
+        }
+        if (err_1 > 0) {
+          y += y_inc;
+          err_1 -= d;
+        }
+        if (err_2 > 0) {
+          z += z_inc;
+          err_2 -= d;
+        }
+        err_0 += dx2;
+        err_1 += dy2;
+        err_2 += dz2;
+        if (err_0 == 0)x += x_inc;
+        if (err_1 == 0)y += y_inc;
+        if (err_2 == 0)z += z_inc;
+    }
+}
+
+void buildMatrix(const Aurora::CudaMatrix &senderList,
+                                  const Aurora::CudaMatrix &receiverList,
+                                  Aurora::CudaMatrix& res,
+                                  const Aurora::CudaMatrix &dims, bool BENT,
+                                  const Aurora::CudaMatrix &potentialMap){
+    unsigned int numDims = senderList.getDimSize(0);
+    unsigned int nTotalRays = senderList.getDimSize(1);
+    unsigned int cnt = 0;
+    unsigned int coeffIndex = 0;
+    unsigned int rayCount = 0; 
+    b3dStruct* structList= nullptr;
+    cudaMalloc((void**)&structList, sizeof(b3dStruct)*nTotalRays);
+    int blocksPerGrid = (nTotalRays + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK;
+    calcLenWeightKernel<<<blocksPerGrid, THREADS_PER_BLOCK>>>(
+        senderList.getData(), receiverList.getData(), res.getData(),structList);
+    unsigned int* offset= nullptr;
+    cudaMalloc((void**)&offset, sizeof(unsigned int)*nTotalRays);
+    calcPathOffsetKernel<<<1,1>>>(structList,offset,nTotalRays);
+    unsigned int size;
+    cudaMemcpy(&size, offset+(nTotalRays-1), 1, cudaMemcpyDeviceToHost);
+    float* iData = nullptr;
+    float* jData = nullptr;
+    float* sData = nullptr;
+
+    cudaMalloc((void**)&iData, sizeof(float)*size);
+    cudaMalloc((void**)&jData, sizeof(float)*size);
+    cudaMalloc((void**)&sData, sizeof(float)*size);
+
+    calcRayBLPath(structList, dims.getData(), offset, jData);
+    cudaFree(structList);
+
+}
+
+
+
--- a/src/transmissionReconstruction/reconstruction/reconstruction.cpp
+++ b/src/transmissionReconstruction/reconstruction/reconstruction.cpp
@@ -8,6 +8,7 @@

 #include "CudaEnvInit.h"
 #include "log/notify.h"
+#include "log/log.h"
 #include "transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h"
 #include "src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.h"

@@ -187,9 +188,12 @@ namespace Recon {

        transParams::nonNeg = false;
        BuildMatrixResult buildMatrixR;
+        
        for(int iter=1; iter<=numIter; ++iter)
        {
+            RECON_INFO("start buildMatrix");
            buildMatrixR = buildMatrix(senderList, receiverList, res, dims, bentRecon && (iter!=1), potentialMap);
+            RECON_INFO("end buildMatrix");
            if(!data.isNull() && bentRecon && iter != numIter)
            {
                //与默认配置bentRecon不符，暂不实现todo
@@ -227,14 +231,19 @@ namespace Recon {
                if (i ==0){
                    if(!data.isNull())
                    {
+                        RECON_INFO("start solveParameterIterator");
                        Matrix sosValue = solveParameterIterator(buildMatrixR.M, b, dims, false, transParams::nonNeg)[0][0];
+                        RECON_INFO("end solveParameterIterator");
                        result.outSOS = slownessToSOS(sosValue, SOS_IN_WATER) ;
+                        RECON_INFO("end slownessToSOS");
                    }
                }
                else{
                    if(!dataAtt.isNull())
                    {
+                        RECON_INFO("start solveParameterIterator");
                        Matrix attValue = solveParameterIterator(buildMatrixR.M, bAtt, dims, false, transParams::nonNeg,1)[0][0];
+                        RECON_INFO("end solveParameterIterator");
                        result.outATT = attValue/100 ;
                    }
                }
--- a/test/Detection_Test.cpp
+++ b/test/Detection_Test.cpp
@@ -8,6 +8,8 @@
 #include "config/config.h"
 #include "common/getMeasurementMetaData.h"
 #include "transmissionReconstruction/detection/detection.h"
+#include "transmissionReconstruction/detection/detection.cuh"
+