Add cuda build matrix

CMakeLists.txt

@@ -9,7 +9,14 @@ 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 "${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)
+endif()
 find_package(Aurora REQUIRED)
 find_package(Parser REQUIRED)
 find_package(Req REQUIRED)
@@ -18,10 +25,25 @@ 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} ${cxx_files} ${header_files} ${Aurora_Source})
-target_compile_options(UR PUBLIC ${Aurora_Complie_Options} "-march=native")
+
+file(GLOB_RECURSE  header_files ./src/*.h ./src/*.cuh)
+
+if (${BUILD_SHARED_LIBS})
+set(cxx_files ./src/UR.cxx)
+
+find_package(MKL CONFIG REQUIRED)
+set(MKL_INTERFACE_FULL intel_lp64)
+add_library(UR SHARED  ${cpp_files} ${cxx_files} ${header_files} ${Aurora_Source})
+target_compile_options(UR PUBLIC $<TARGET_PROPERTY:MKL::MKL,INTERFACE_COMPILE_OPTIONS>)
+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)
+else()
+set(cxx_files ./src/main.cxx)
+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/)
 target_include_directories(UR PUBLIC ${Aurora_INCLUDE_DIRS})
 target_include_directories(UR PUBLIC ${DCMTK_INCLUDE_DIRS})
@@ -51,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)
 
@@ -62,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)

src/UR.cxx

@@ -0,0 +1,67 @@
+#include "Parser.h"
+#include "Parser.h"
+#include "config/config.h"
+#include "log/notify.h"
+#include "log/notify.h"
+#include <exception>
+#include <iostream>
+#include <vector>
+#include <stdexcept>
+#define EIGEN_USE_MKL_ALL
+
+#include "src/common/fileHelper.h"
+#include "startReconstructions.h"
+#include "log/log.h"
+extern "C" {
+    int ReconProcess(const char* reconID,const char* path,const char* refPath,const char* output,const char* config){
+        try{
+        int argNum = 5; 
+        std::vector<std::string> args(argNum);
+        std::string outPutPath = output?output:Recon::DEFAULT_OUTPUT_PATH;
+        outPutPath = Recon::fixPathSlash(outPutPath);
+        std::string directoryPath = outPutPath;
+        if(!Recon::isDirectory(directoryPath))
+        {
+            printf("Output directory is not valid.");
+            return -4;
+        }
+        auto defaultLogger = getLogger("Main",outPutPath.data());
+        spdlog::set_default_logger(defaultLogger);
+        if(!reconID)
+        {
+            RECON_INFO("No recon id.");
+            return -1;
+        }
+        std::string ReconID = reconID;
+        if(!path)
+        {
+            RECON_INFO("No reconstruction data.");
+            return -2;
+        }
+        std::string configPath = Recon::fixPathSlash(config?config:Recon::DEFAULT_CONFIG_PATH);
+        Recon::initalizeConfig(configPath);
+        if( !refPath && Recon::transParams::runTransmissionReco)
+        {
+            RECON_INFO("Running transmission reconstruction, but no refrence data.");
+            return -3;
+        }
+        std::string dataPath = Recon::fixPathSlash(path);
+        std::string dataRefPath = Recon::fixPathSlash(refPath);
+        RECON_INFO("start");
+        Recon::notifyStart("");
+        int exitcode = Recon::startReconstructions(dataPath, dataRefPath, outPutPath);
+        if (exitcode == 0)
+        {
+            RECON_INFO("finish");
+            return exitcode;
+        }
+        else{
+            return exitcode;
+        }
+        }
+        catch(const std::exception& ex){
+            RECON_INFO("Recon fail by unknow{}", ex.what());
+            return -99;
+        }
+    }
+}

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:
             {

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);

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;

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());

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;
         }
     }
 }

src/transmissionReconstruction/detection/detection.cu

@@ -0,0 +1,326 @@
+#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));
+    cudaDeviceSynchronize();
+    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);
+    cudaDeviceSynchronize();
+    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));
+        }
+        cudaDeviceSynchronize();
+        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);
+        cudaDeviceSynchronize();
+        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;
+    }

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__

src/transmissionReconstruction/reconstruction/buildMatrix/Bresenham.cpp

@@ -1,4 +1,5 @@
 #include "Bresenham.h"
+#include <cstdio>
 #include <math.h> 
 #include <stdlib.h>
 #include <vector>
@@ -109,13 +110,89 @@ 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;
+    int f = -1;
+    if ((l >= m) && (l >= n)){
+        maxV = l;
+        d = dx2;
+        dx2 = 0;
+        f =0;
+    }
+    else if((m >= l) && (m >= n)){
+        maxV = m;
+        d = dy2;
+        dy2 = 0;
+        f = 1;
+    }
+    else{
+        maxV = n;
+        d = dz2;
+        dz2 = 0;
+        f =2;
+    }
+    err_0 = f==0?0:(dx2 - maxV);
+    err_1 = f==1?0:(dy2 - maxV);
+    err_2 = f==2?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 (f == 0)x += x_inc;
+        if (f == 1)y += y_inc;
+        if (f == 2)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;
     auto temp  =
         Recon::b3dMexDouble(StartPt.getData(), endPts.getData(), size);
-    return Aurora::Matrix::fromRawData(temp, 3, size / 3);
+        return Aurora::Matrix::fromRawData(temp, 3, size / 3);
 }
 
 float * b3dMexDouble(float *startPoint,float *endPoint, size_t& outputSize) {
@@ -215,5 +292,7 @@ float * b3dMexDouble(float *startPoint,float *endPoint, size_t& outputSize) {
     outputSize = points.size();
     return output;
 }
+
+
 }
 

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);

src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cu

@@ -0,0 +1,238 @@
+#include "Function2D.cuh"
+
+
+#include <cstdio>
+#include <thrust/reduce.h>
+#include <thrust/execution_policy.h>
+#include <cuda_runtime.h>
+#include <math.h>
+
+#include "CudaMatrix.h"
+#include "buildMatrix.cuh"
+#include "Sparse.h"
+
+
+namespace  {
+    const int THREADS_PER_BLOCK = 256;
+}
+
+
+
+__device__ float getPixelLengthKernel(float* aStartPt, float*aEndPt, float* aRes, int aPathLenDisc)
+{
+    float temp1;
+    float temp2;
+    float temp3;
+
+    temp1 = aEndPt[0]- aStartPt[0];
+    temp2 = aEndPt[1]- aStartPt[1];
+    temp3 = aEndPt[2]- aStartPt[2];
+    temp1 *= aRes[0];
+    temp2 *= aRes[1];
+    temp3 *= aRes[2];
+    temp1 = temp1*temp1;
+    temp2 = temp2*temp2;
+    temp3 = temp3*temp3;
+    return (sqrtf(temp3+temp2+temp1))/(float)aPathLenDisc;
+}
+
+__global__ struct b3dStruct{
+    int max, f, pathlen;
+    float  x_inc, y_inc, z_inc,err_0, err_1, err_2, dx2, dy2, dz2, d, x, y, z, weight;
+};
+__device__ void getRayBLen(float *startPoint, float *endPoint, 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;
+    float vmax ;
+    if ((l >= m) && (l >= n)){
+        aOut.max =ceilf(l);
+        vmax = l;
+        aOut.d = aOut.dx2;
+        aOut.dx2 = 0;
+        aOut.f = 0;
+    }
+    else if((m >= l) && (m >= n)){
+        aOut.max = ceilf(m);
+        vmax = m;
+        aOut.d = aOut.dy2;
+        aOut.dy2 = 0;
+        aOut.f = 1;
+    }
+    else{
+        aOut.max = ceilf(n);
+        vmax = n;
+        aOut.d = aOut.dz2;
+        aOut.dz2 = 0;
+        aOut.f = 2;
+    }
+    aOut.err_0 = aOut.f==0?0:(aOut.dx2 - vmax);
+    aOut.err_1 = aOut.f==1?0:(aOut.dy2 - vmax);
+    aOut.err_2 = aOut.f==2?0:(aOut.dz2 - vmax);
+}
+
+/**
+ * @brief 计算traceStraightRayBresenham的pathlen和weight的核函数
+ * 
+ * @param aStratPt 
+ * @param aEndPt 
+ * @param aSize 
+ * @param aRes 
+ * @param aOuts 
+ * @return _
+ */
+__global__ void calcLenWeightKernel(float* aStratPt, float* aEndPt, float* aRes,int aSize, 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,aOuts[idx].max);
+};
+
+__global__ void calcPathOffsetKernel(b3dStruct* aStructs ,unsigned int* aPath, int aSize){
+    aPath[0] = aStructs[0].max;
+    for (int i=1; i<aSize; i++) {
+        aPath[i] =aStructs[i].max+aPath[i-1]; 
+    }
+}
+
+__device__ float convertToLinearIndices(const float3& aDims, const float3& aPath){
+    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;
+    return coord_col1 + coord_col2Sub1 * matrixSize_1 + coord_col3Sub1 * matrixSize_1 * matrixSize_2 - 1 ;
+}
+
+__global__ void calcRayBLPathKernel(struct b3dStruct* aInstruct,
+            float* dims, unsigned int * aOffset,
+            float* aJOut, unsigned int aSize) { 
+    unsigned int idx = blockIdx.x;
+    if (idx > aSize)return;
+    b3dStruct* in = aInstruct+idx;
+    int f;
+    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;
+
+    dx2 = in->dx2;
+    dy2 = in->dy2;
+    dz2 = in->dz2;
+    d = in->d;
+
+    err_0 = in->err_0;
+    err_1 = in->err_1;
+    err_2 = in->err_2;
+    max = in->max;
+    f = in->f;
+
+    float3 dims3{dims[0],dims[1],dims[2]};
+    float* output_ptr = aJOut + (idx ==0?0:aOffset[idx-1]);
+    for (float i = 0; i < max; i++) {
+        output_ptr[0]=convertToLinearIndices(dims3, {roundf(x),roundf(y),roundf(z)});
+        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 (f == 0) x += x_inc;
+        if (f == 1) y += y_inc;
+        if (f == 2) z += z_inc;
+    }
+}
+
+__global__ void fillIAndSKernel(struct b3dStruct* aInstructs,unsigned int * aOffset,float* aIOut,float* aSOut){
+    b3dStruct& s = aInstructs[blockIdx.x];
+    __shared__ unsigned int beginIdx; 
+    __shared__ unsigned int length;
+    __shared__ float weight; 
+
+    if (threadIdx.x == 0){
+        beginIdx = (blockIdx.x == 0) ? 0 : aOffset[blockIdx.x - 1];
+        length = s.max;
+        weight = s.weight;
+    }
+    __syncthreads();
+
+    for (int i = 0; threadIdx.x + i < length; i+=blockDim.x) {
+        aIOut[beginIdx + threadIdx.x + i] = blockIdx.x; 
+        aSOut[beginIdx + threadIdx.x + i] = weight; 
+    }
+}
+
+Recon::BuildMatrixResult Recon::buildMatrix(const Aurora::CudaMatrix &senderList,
+                                  const Aurora::CudaMatrix &receiverList,
+                                  Aurora::CudaMatrix& res,
+                                  const Aurora::CudaMatrix &dims, bool BENT,
+                                  const Aurora::CudaMatrix &potentialMap){
+    Recon::BuildMatrixResult result;
+    unsigned int numDims = senderList.getDimSize(0);
+    unsigned int nTotalRays = senderList.getDimSize(1);
+    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(),nTotalRays,structList);
+    cudaDeviceSynchronize();
+
+    unsigned int* offset= nullptr;
+    cudaMalloc((void**)&offset, sizeof(unsigned int)*nTotalRays);
+    calcPathOffsetKernel<<<1,1>>>(structList,offset,nTotalRays);
+    cudaDeviceSynchronize();
+
+    unsigned int size;
+    cudaMemcpy(&size, offset+(nTotalRays-1), sizeof(unsigned int), 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);
+
+    calcRayBLPathKernel<<<nTotalRays,1>>>(structList, dims.getData(), offset, jData, nTotalRays);
+    fillIAndSKernel<<<nTotalRays,THREADS_PER_BLOCK>>>(structList, offset, iData, sData);
+    cudaDeviceSynchronize();
+    cudaFree(structList);
+
+    auto mI = Aurora::CudaMatrix::fromRawData(iData, 1, size ,1);
+    auto mJ = Aurora::CudaMatrix::fromRawData(jData, 1, size,1);
+    auto mS = Aurora::CudaMatrix::fromRawData(sData, 1, size,1);
+    result.M = Aurora::Sparse(mI.toHostMatrix(), mJ.toHostMatrix(),mS.toHostMatrix(),nTotalRays, Aurora::prod(dims).getValue(0));
+    return result;   
+
+}
+
+
+

src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cuh

@@ -0,0 +1,12 @@
+#ifndef __BUILDMATRIX_CUH__
+#define __BUILDMATRIX_CUH__
+#include "buildMatrix.h"
+namespace Recon {
+    
+BuildMatrixResult buildMatrix(const Aurora::CudaMatrix &senderList,
+                                  const Aurora::CudaMatrix &receiverList,
+                                  Aurora::CudaMatrix& res,
+                                  const Aurora::CudaMatrix &dims, bool BENT,
+                                  const Aurora::CudaMatrix &potentialMap);
+}
+#endif // __BUILDMATRIX_H__

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
@@ -222,19 +226,24 @@ namespace Recon {
             {
                 allHitMaps.push_back(buildMatrixR.hitmap);
             }
-            #pragma omp parallel for num_threads(2)
+            // #pragma omp parallel for num_threads(2)
             for (int i =0; i<2; i++){
                 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 ;
                     }
                 }

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"
+