feat: modify getTransmissionData function , use detectTofThreshold function to do get TransmissionData

src/transmissionReconstruction/detection/detection.cpp

@@ -82,11 +82,12 @@ namespace Recon {
     }
     
     TimeWindowResult applyTimeWindowing(const Aurora::Matrix &AscanBlock, float sampleRate,
-        const Aurora::Matrix &distBlock, float sosWater,
+        const Aurora::Matrix &distBlock, const Aurora::Matrix& sosBlock,float sosWater,
         float startOffset, float segmentLenOffset,
         float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow)
     {
-        Aurora::Matrix sosOffset = Aurora::Matrix::fromRawData(new float[1]{0}, 1);
+        // Aurora::Matrix sosOffset = Aurora::Matrix::fromRawData(new float[1]{0}, 1);
+        Aurora::Matrix sosOffset = calculateSOSOffset(sosBlock, sosWater, distBlock, sampleRate);
         auto calcResult = calculateStarEndSearchPosition(distBlock, minSpeedOfSound, maxSpeedOfSound, sampleRate,AscanBlock.getDimSize(0), sosOffset, startOffset, segmentLenOffset);
 
         auto AscanBlockProcessed = zeros(AscanBlock.getDimSize(0),AscanBlock.getDimSize(1));
@@ -173,6 +174,7 @@ namespace Recon {
     DetectResult detectTofVectorized(
         const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
         const Aurora::Matrix &distBlock, const Aurora::Matrix &distBlockRef,
+        const Aurora::Matrix &sosWaterBlock,const Aurora::Matrix &sosWaterRefBlock,
         float aSOSWater, float aSOSWaterRef,
         int useTimeWindowing, int aScanReconstructionFrequency,
         float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound,
@@ -187,11 +189,11 @@ namespace Recon {
         timeResult2.AscanBlockProcessed = AscanRefBlock;
         if (useTimeWindowing == 1) {
             timeResult1 = applyTimeWindowing(
-                AscanBlock, sampleRate, distBlock, 
+                AscanBlock, sampleRate, distBlock, sosWaterBlock,
                 aSOSWater, offsetElectronicSamples, detectionWindowSOS,
                 minSpeedOfSound, maxSpeedOfSound, gaussWindow);
             timeResult2 = applyTimeWindowing(
-                AscanRefBlock, sampleRate, distBlockRef, 
+                AscanRefBlock, sampleRate, distBlockRef, sosWaterRefBlock,
                 aSOSWaterRef, offsetElectronicSamples, detectionWindowSOS,
                 minSpeedOfSound, maxSpeedOfSound, gaussWindow);
 
@@ -244,6 +246,7 @@ namespace Recon {
     DetectResult detectTofAndAtt(
         const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
         const Aurora::Matrix &distBlock, const Aurora::Matrix &distRefBlock,
+        const Aurora::Matrix &sosWaterBlock,const Aurora::Matrix &sosWaterRefBlock,
         int resampleFactor,int nthreads, float aSOSWater, float aSOSWaterRef,
         int useTimeWindowing, int aScanReconstructionFrequency,int detectionWindowATT,
         float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound,
@@ -258,11 +261,11 @@ namespace Recon {
         timeResult2.AscanBlockProcessed = AscanRefBlock;
         if (useTimeWindowing == 1) {
             timeResult1 = applyTimeWindowing(
-                AscanBlock, sampleRate, distBlock, 
+                AscanBlock, sampleRate, distBlock, sosWaterBlock,
                 aSOSWater, offsetElectronicSamples, detectionWindowSOS,
                 minSpeedOfSound, maxSpeedOfSound, gaussWindow);
             timeResult2 = applyTimeWindowing(
-                AscanRefBlock, sampleRate, distRefBlock, 
+                AscanRefBlock, sampleRate, distRefBlock, sosWaterRefBlock,
                 aSOSWaterRef, offsetElectronicSamples, detectionWindowSOS,
                 minSpeedOfSound, maxSpeedOfSound, gaussWindow);
 
@@ -420,7 +423,7 @@ namespace Recon {
                                        const Aurora::CudaMatrix &AscanRefBlock,
                                        const Aurora::CudaMatrix &distBlock,
                                        const Aurora::CudaMatrix &distRefBlock,
-                                       float aSOSWater, float aSOSWaterRef)
+                                       float aSOSWater, float aSOSWaterRef, float aOffset)
     {
         switch (Recon::transParams::version) {
         // case 1: {
@@ -433,19 +436,40 @@ namespace Recon {
         //         Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow);
         // }
         // case 2:
+        // {
+        //     auto SOSWaterBlock = zerosCuda(1,AscanBlock.getDimSize(1));
+        //     SOSWaterBlock.setBlockValue(1,0,AscanBlock.getDimSize(1)-1,aSOSWater);
+        //     auto SOSWaterRefBlock = zerosCuda(1,AscanBlock.getDimSize(1));
+        //     SOSWaterRefBlock.setBlockValue(1,0,AscanBlock.getDimSize(1)-1,aSOSWaterRef);
+        //     auto r =  detectTofAndAtt(
+        //         AscanBlock, AscanRefBlock, distBlock, distRefBlock,SOSWaterBlock,SOSWaterRefBlock,
+        //         Recon::transParams::resampleFactor, Recon::transParams::nThreads,
+        //         aSOSWater, aSOSWaterRef, 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;
+        // }
+        // case 3:
         default:
-            auto r =  detectTofAndAtt(
+        {
-                AscanBlock, AscanRefBlock, distBlock, distRefBlock,
+            float sampleRate = Recon::transParams::aScanReconstructionFrequency;
-                Recon::transParams::resampleFactor, Recon::transParams::nThreads,
+            auto tof = detectTofThreshold(AscanBlock, distBlock, aOffset, aSOSWater, sampleRate, 
-                aSOSWater, aSOSWaterRef, Recon::transParams::useTimeWindowing,
+                Recon::transParams::offsetElectronic*sampleRate, 
-                Recon::transParams::aScanReconstructionFrequency, Recon::transParams::detectionWindowATT,
+                Recon::transParams::gaussWindow, Recon::transParams::minSpeedOfSound,
-                Recon::transParams::offsetElectronic, Recon::transParams::detectionWindowSOS, Recon::transParams::minSpeedOfSound,
+                Recon::transParams::maxSpeedOfSound);
-                Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow);
+            auto att = detectAttVectorizedCuda(AscanBlock, AscanRefBlock,distRefBlock, aSOSWaterRef, tof,
+                    sampleRate, Recon::transParams::offsetElectronic,
+                    Recon::transParams::detectionWindowATT);
             DetectResult ret;
-            ret.att = r.att.toHostMatrix();
+            ret.tof = (tof - distBlock/aSOSWater).toHostMatrix();
-            ret.sosValue = r.sosValue.toHostMatrix();
+            ret.att = att.toHostMatrix();
-            ret.tof = r.tof.toHostMatrix();
             return ret;
         }
+        }
     }
 }

src/transmissionReconstruction/detection/detection.cu

@@ -99,6 +99,115 @@ CudaMatrix Recon::calculateSOSOffset(const CudaMatrix &sosBlock, float reference
     return offset;
 }
 
+__global__ void addWinKernel(float* aData, int aRowCount, int aColCount, float aOffset, float aMinSOS, float aMaxSOS, float* aDist, float aSamepleRate)
+{
+    __shared__ int startIndex;
+    __shared__ int endIndex;
+    unsigned int i = blockIdx.x;
+    unsigned int base_offset = i * aRowCount;
+    float * dataPointer = aData + base_offset;
+    float dist = aDist[i];
+    //确定每个thread的index
+    if (threadIdx.x == 0)
+    {
+        float offset = aOffset * aSamepleRate;
+        startIndex = (int)(round(dist / aMaxSOS * aSamepleRate) + offset)-1;
+        endIndex = (int)(round(dist / aMinSOS * aSamepleRate) + offset)-1;
+    }
+    for (size_t j = threadIdx.x; j < aRowCount; j+=blockDim.x)
+    {
+        if(j<startIndex) dataPointer[j] = .0;
+        if(j>endIndex) dataPointer[j] = .0;
+    }
+}
+
+__global__ void genTofKernel(float* aData, int aRowCount, int aColCount, float aOffset, float* aTOFRef, float aSamepleRate, int* index, int aNpeaks)
+{
+    unsigned int i = blockIdx.x;
+    unsigned int base_offset = i * aRowCount ;
+    float * dataPointer = aData + base_offset;
+    int startIndex = index[i*aNpeaks]-15;
+    startIndex = startIndex < 0 ? 0 : startIndex;
+    int endIndex = index[i*aNpeaks]+15;
+    endIndex = endIndex >= aRowCount? aRowCount-1 : endIndex;
+    int maxIndex = 0;
+    float maxValue = FLT_MIN;
+    //find max in index[i] + (-15,15)
+    for (size_t j = startIndex; j <= endIndex; j++)
+    {
+        if (dataPointer[j] > maxValue){
+            maxIndex = j;
+            maxValue = dataPointer[j];
+        }
+    }
+    index[i] = maxIndex;
+    //attention: this remove redudent -aTOFRef[i] ,because it will be add back in upper called level
+    aTOFRef[i] = ((float)(maxIndex+1))/aSamepleRate + aOffset ;
+    // aTOFRef[i] = ((float)(maxIndex+1))/aSamepleRate + aOffset - aTOFRef[i];
+
+}
+
+CudaMatrix Recon::detectTofThreshold(const CudaMatrix& aAscan, const CudaMatrix &aDistRef, float aOffset, float aSOSWater,
+            float aAScanReconstructionFrequency, float aOffsetElectronic, int aUseTimeWindowing, float aMinSpeedOfSound,
+            float aMaxSpeedOfSound)
+{
+    const int THREADS_PER_BLOCK_ = 32;
+    float offsetWV = -aOffset;
+    float sampleRate = aAScanReconstructionFrequency;
+    float offsetElectronSamples = aOffsetElectronic * sampleRate;
+    // if (inits.detection.useTimeWindowing == 1)
+	// parfor idx__ = 1 : size(AscanBlock, 2)
+    //     AscanBlock(:, idx__)     	= add_win(AscanBlock(:, idx__) ,...
+    //                                           offsetElectronicSamples,...
+    //                                           inits.detection.minSpeedOfSound,...
+    //                                           inits.detection.maxSpeedOfSound,...
+    //                                           distBlock(idx__), sampleRate);
+	//     end
+    // end
+    int * index = nullptr;
+    if (aUseTimeWindowing){
+        addWinKernel<<<aAscan.getDimSize(1), THREADS_PER_BLOCK_>>>(aAscan.getData(), aAscan.getDimSize(0),
+                aAscan.getDimSize(1), offsetElectronSamples, 
+                aMinSpeedOfSound, aMaxSpeedOfSound, aDistRef.getData(), sampleRate);
+        cudaDeviceSynchronize();
+        
+    }
+    // AScan_env                                   = envelope(AscanBlock);
+    // pkValue_mea                                 = max(AScan_env, [], 1);
+    // AScan_env_norm                              = AScan_env ./ pkValue_mea;
+    // parfor idx__ = 1 : size(AScan_env_norm, 2)
+    //     %-find 1st wvfrom based on envelope
+    //     [~, locs]                               = findpeaks(AScan_env_norm(:, idx__),...
+    //                                                         'MinPeakHeight', 0.2,...
+    //                                                         'Npeaks', 10,...
+    //                                                         'MInPeakProminence', 0.05);
+    //     tof_env                                 = locs(1);
+    //     %-then select maximum pt on match-filted signal
+    //     Ascan_valid_part                        = zeros(size(AscanBlock, 1), 1);
+    //     Ascan_valid_part((-15:15) + tof_env)    = AscanBlock((-15:15) + tof_env, idx__);
+    //     loc                                     = find(Ascan_valid_part == max(Ascan_valid_part), 1);
+    //     try
+    //     tof_mea(idx__)                          = loc;
+    //     catch ME
+    //         error(ME)
+    //     end
+    // end
+    {
+        auto envAScanTemp = abs(hilbert(aAscan));//envelope(Ascanblock);
+        auto pkValueMea = max(envAScanTemp);
+        envAScanTemp = envAScanTemp / pkValueMea;
+        auto pearks = Aurora::findPeaks(envAScanTemp, 10, 0.2, 0.05, &index);
+    }
+
+    auto tofRef =  aDistRef.deepCopy();
+    genTofKernel<<<aAscan.getDimSize(1),1>>>(aAscan.getData(), aAscan.getDimSize(0),
+            aAscan.getDimSize(1), offsetWV, tofRef.getData(),
+            sampleRate,index, 10);
+    cudaDeviceSynchronize();
+    cudaFree(index);
+    return tofRef;
+}
+
 Recon::SearchPositionC Recon::calculateStarEndSearchPosition(const CudaMatrix &aVDistBlock,
                                 float minSpeedOfSound, float maxSpeedOfSound,
                                 float sampleRate, float maxSample,
@@ -170,11 +279,13 @@ __global__ void guassWindowKernel(float* aStartSearch,float* aEndSearch,
 }
 
 Recon::TimeWindowResultC Recon::applyTimeWindowing(const Aurora::CudaMatrix &AscanBlock, float sampleRate,
-        const Aurora::CudaMatrix &distBlock, 
+        const Aurora::CudaMatrix &distBlock,const Aurora::CudaMatrix &sosWaterBlock,
         float aSOSWater, float startOffset, float segmentLenOffset,
         float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow)
     {
-        Aurora::CudaMatrix sosOffset = Aurora::zerosCuda(1,1);
+        // Aurora::CudaMatrix sosOffset = Aurora::zerosCuda(1,1);
+        Aurora::CudaMatrix sosOffset = calculateSOSOffset(sosWaterBlock,aSOSWater, distBlock, sampleRate);
+
         auto calcResult = calculateStarEndSearchPosition(distBlock, minSpeedOfSound, maxSpeedOfSound, sampleRate,AscanBlock.getDimSize(0), sosOffset, startOffset, segmentLenOffset);
 
         auto AscanBlockProcessed = zerosCuda(AscanBlock.getDimSize(0),AscanBlock.getDimSize(1));
@@ -253,6 +364,7 @@ int nextpow2(unsigned int value){
 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 aSOSWater, float aSOSWaterRef,
         int useTimeWindowing, int aScanReconstructionFrequency,int detectionWindowATT,
         float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound,
@@ -267,11 +379,11 @@ Recon::DetectResultC Recon::detectTofAndAtt(
         timeResult2.AscanBlockProcessed = AscanRefBlock;
         if (useTimeWindowing == 1) {
             timeResult1 = applyTimeWindowing(
-                AscanBlock, sampleRate, distBlock, 
+                AscanBlock, sampleRate, distBlock, soswaterBlock,
                 aSOSWater, offsetElectronicSamples, detectionWindowSOS,
                 minSpeedOfSound, maxSpeedOfSound, gaussWindow);
             timeResult2 = applyTimeWindowing(
-                AscanRefBlock, sampleRate, distRefBlock, 
+                AscanRefBlock, sampleRate, distRefBlock, soswaterRefBlock,
                 aSOSWaterRef, offsetElectronicSamples, detectionWindowSOS,
                 minSpeedOfSound, maxSpeedOfSound, gaussWindow);
 

src/transmissionReconstruction/detection/detection.cuh

@@ -30,7 +30,7 @@ struct TimeWindowResultC {
                         float offsetElectronic, int detectionWindowATT);
     
     TimeWindowResultC applyTimeWindowing(const Aurora::CudaMatrix &AscanBlock, float sampleRate,
-        const Aurora::CudaMatrix &distBlock,
+        const Aurora::CudaMatrix &distBlock,const Aurora::CudaMatrix &sosWaterBlock,
         float aSOSWater, float startOffset, float segmentLenOffset,
         float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow);
     SearchPositionC calculateStarEndSearchPosition(const CudaMatrix &aVDistBlock,
@@ -39,9 +39,14 @@ struct TimeWindowResultC {
                                 const CudaMatrix &aVSosOffsetBlock,
                                 float startOffset, float segmentLenOffset);
     CudaMatrix calculateSOSOffset(const CudaMatrix &sosBlock, float referenceSOS, const CudaMatrix &distBlock, float sampleRate);
+
+    CudaMatrix detectTofThreshold(const CudaMatrix& aAscan, const CudaMatrix &aDistRef,float aOffset, float aSOSWater,
+            float aAScanReconstructionFrequency, float aOffsetElectronic, int aUseTimeWindowing, float aMinSpeedOfSound, float aMaxSpeedOfSound);
+
     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 aSOSWater, float aSOSWaterRef,
         int useTimeWindowing, int aScanReconstructionFrequency,int detectionWindowATT,
         float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound,

src/transmissionReconstruction/detection/detection.h

@@ -32,8 +32,8 @@ calculateStarEndSearchPosition(const Aurora::Matrix &aVDistBlock,
 
 TimeWindowResult applyTimeWindowing(
     const Aurora::Matrix &AscanBlock, float sampleRate,
-    const Aurora::Matrix &distBlock, float sosWater,
+    const Aurora::Matrix &distBlock, const Aurora::Matrix &sosWaterBlock,
-    float startOffset, float segmentLenOffset,
+    float sosWater, float startOffset, float segmentLenOffset,
     float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow);
 
 Aurora::Matrix
@@ -56,6 +56,7 @@ DetectResult
 detectTofAndAtt(
     const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
     const Aurora::Matrix &distBlock, const Aurora::Matrix &distRefBlock,
+    const Aurora::Matrix &sosWaterBlock,const Aurora::Matrix &sosWaterRefBlock,
     int resampleFactor, int nthreads, float aSOSWater, float aSOSWaterRef,
     int useTimeWindowing, int aScanReconstructionFrequency,
     int detectionWindowATT, float offsetElectronic, int detectionWindowSOS,
@@ -77,7 +78,7 @@ DetectResult
 transmissionDetection(
     const Aurora::CudaMatrix &AscanBlock, const Aurora::CudaMatrix &AscanRefBlock,
     const Aurora::CudaMatrix &distBlock, const Aurora::CudaMatrix &distRefBlock,
-    float aSOSWater, float aSOSWaterRef);
+    float aSOSWater, float aSOSWaterRef, float aOffset);
 
 } // namespace Recon
 

src/transmissionReconstruction/detection/getTransmissionData.cpp

@@ -288,7 +288,7 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
                 cudaSetDevice(index % BUFFER_SIZE); 
                 DetectResult detect = transmissionDetection( blockData.ascanBlock, blockData.ascanBlockRef, 
                                                              blockData.dists.toDeviceMatrix(), blockData.distRefBlock.toDeviceMatrix(), 
-                                                             aSOS.expectedSOSWater[0], aSOSRef.expectedSOSWater[0]);
+                                                             aSOS.expectedSOSWater[0], aSOSRef.expectedSOSWater[0], aPreComputes.offset);
                 blockData.attData = detect.att;
                 blockData.tofData = detect.tof;
                 BlockOfTransmissionData transmissionBlock=blockData;