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

src/common/SoundSpeedWaterCalculation/estimateSOSWater.cpp

@@ -0,0 +1,143 @@
+#include "estimateSOSWater.h"
+#include "Function1D.cuh"
+#include "Function1D.h"
+#include "Function3D.h"
+#include "Matrix.h"
+#include "CudaMatrix.h"
+#include "common/dataBlockCreation/getAscanBlock.h"
+#include "common/getMeasurementMetaData.h"
+#include "transmissionReconstruction/detection/getTransmissionData.cuh"
+#include "transmissionReconstruction/detection/detection.h"
+#include "config/config.h"
+#include "config/geometryConfig.h"
+
+#include "Function2D.cuh"
+#include "Function2D.h"
+#include "common/dataBlockCreation/getAscanBlock.h"
+#include "common/dataBlockCreation/removeDataFromArrays.h"
+#include "log/log.h"
+#include <cstring>
+
+using namespace Recon;
+
+namespace
+{
+    const Aurora::Matrix EMIT_TAS = Aurora::Matrix::fromRawData(new float[17]{94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110}, 17);
+    const Aurora::Matrix RECEIVE_TAS = Aurora::Matrix::fromRawData(new float[17]{103, 104, 105, 106, 107, 108, 109, 110, 94, 95, 96, 97, 98, 99, 100, 101, 102}, 17);
+    const Aurora::Matrix ALL_RECEIVER_TAS_LIST = Aurora::Matrix::fromRawData(new float[128] {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128},128); //1~128
+    const Aurora::Matrix ALL_RECEIVER_ELEMENT_LIST = Aurora::Matrix::fromRawData(new float[18] {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18},18);
+    const float THRE = 0.5;
+    inline int findIndexFromEmitterAndReceiver(const Aurora::Matrix& aEmitter, float aEmitterValue, const Aurora::Matrix& aReceiver, float aReceiverValue)
+    {
+        for (int i = 0; i < aEmitter.getDataSize(); ++i)
+        {
+            if (aEmitter[i] == aEmitterValue && aReceiver[i] == aReceiverValue)
+            {
+                return i;
+            }
+        }
+        return -1;
+    }
+
+    Aurora::Matrix matchFilt(const Aurora::CudaMatrix& aSig, const Aurora::CudaMatrix& aMod)
+    {
+        Aurora::CudaMatrix fftMod = fft(aMod);
+        Aurora::CudaMatrix fftSig = fft(aSig);
+        Aurora::CudaMatrix mfFftSig = getTransmissionDataSubFunction(fftSig, fftMod);
+        Aurora::Matrix mfSig = Aurora::real(ifft(mfFftSig)).toHostMatrix();
+        std::memmove(mfSig.getData(), mfSig.getData() + 1, (mfSig.getDataSize() - 1) * sizeof(float));
+        return mfSig;
+    }
+
+    Aurora::Matrix addWin(const Aurora::Matrix& aSig, float aOffset, float aMinSos, float aMaxSos, float aDist, float aFq)
+    {
+        float offset = aOffset *  aFq;
+        int start_ = round(aDist / aMaxSos * aFq);
+        int end_ = round(aDist / aMinSos * aFq);
+        float start = start_ + offset;
+        float end = end_ + offset;
+        Aurora::Matrix result = Aurora::zeros(aSig.getDimSize(0), aSig.getDimSize(1));
+        result.setBlock(0, start-1, end-1, aSig.block(0, start-1, end-1));
+        return result;
+    }
+
+    inline int findFirstvalueGreaterThanGivenValue(const Aurora::Matrix& aMatrix, float aValue)
+    {
+        for(int i=0; i<aMatrix.getDataSize(); ++i)
+        {
+            if(aMatrix[i] > aValue)
+            {
+                return i;
+            }
+        }
+        return -1;
+    }
+}
+
+Aurora::Matrix Recon::estimateSOSWater(Parser *aParser, const CEInfo &aCE, const PreComputes& aPreComputes)
+{
+    float offset = -aPreComputes.offset;
+    float thre = 0.99;
+    Aurora::Matrix tof = Aurora::zeros(1, EMIT_TAS.getDataSize() * 18 * 18);
+    Aurora::Matrix distance = tof.deepCopy();
+    Aurora::Matrix mp = Aurora::Matrix::fromRawData(new float[1]{1}, 1);
+    Aurora::Matrix snMatrix = Aurora::Matrix::fromRawData(new float[18]{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18}, 18);
+
+    for(int idx=0; idx<EMIT_TAS.getDataSize(); ++idx)
+    {
+        Aurora::Matrix slMatrix = Aurora::Matrix::fromRawData(new float[1]{EMIT_TAS[idx]}, 1);
+        AscanBlock ascanTotal = getAscanBlock(aParser, mp, slMatrix, snMatrix, ALL_RECEIVER_TAS_LIST, ALL_RECEIVER_ELEMENT_LIST);
+        #pragma omp parallel for
+        for(int sn=1; sn<=18; ++sn)
+        {
+            // Aurora::Matrix snMatrix = Aurora::Matrix::fromRawData(new float[1]{sn}, 1);
+            // Aurora::Matrix slMatrix = Aurora::Matrix::fromRawData(new float[1]{EMIT_TAS[idx]}, 1);
+            // AscanBlock ascan = getAscanBlock(aParser, mp, slMatrix, snMatrix, ALL_RECEIVER_TAS_LIST, ALL_RECEIVER_ELEMENT_LIST);
+            // size_t dataSize = ascan.ascanBlock.getDataSize();
+            // short* data = new short[dataSize];
+            // std::copy(ascan.ascanBlock.getData(), ascan.ascanBlock.getData() + dataSize, data);
+            // Aurora::Matrix rfSet = Aurora::convertfp16tofloat(data, ascan.ascanBlock.getDimSize(0), ascan.ascanBlock.getDimSize(1));
+            // delete[] data;
+            // Aurora::CudaMatrix rfSetDevice = rfSet.toDeviceMatrix() / ascan.gainBlock.toDeviceMatrix();
+            //AscanBlock ascan = getAscanBlock(aParser, mp, slMatrix, snMatrix, ALL_RECEIVER_TAS_LIST, ALL_RECEIVER_ELEMENT_LIST);
+            size_t dataSize = ascanTotal.ascanBlock.getDataSize()  / 18;
+            short* data = new short[dataSize];
+            std::copy(ascanTotal.ascanBlock.getData() + (sn - 1) * dataSize, ascanTotal.ascanBlock.getData() + sn * dataSize, data);
+            Aurora::Matrix rfSet = Aurora::convertfp16tofloat(data, ascanTotal.ascanBlock.getDimSize(0), ascanTotal.ascanBlock.getDimSize(1) / 18);
+            delete[] data;
+            Aurora::CudaMatrix rfSetDevice = rfSet.toDeviceMatrix() / ascanTotal.gainBlock.block(1, sn-1, sn-1).toDeviceMatrix();
+            for(int rn=1; rn<=18; ++rn)
+            {
+                int cnt = idx * 18 * 18 + 18 * (sn - 1) + rn - 1;
+                Aurora::Matrix positionOfEmitter = emitterPositions[EMIT_TAS[idx]-1](sn-1).toMatrix();
+                Aurora::Matrix positionOfReceiver = emitterPositions[RECEIVE_TAS[idx]-1](rn-1).toMatrix();
+                distance[cnt] = Aurora::norm(positionOfEmitter - positionOfReceiver, Aurora::Norm2);
+                int rfIndex = findIndexFromEmitterAndReceiver(aCE.tasIndices, RECEIVE_TAS[idx], aCE.receiverIndices, rn);
+                Aurora::CudaMatrix rf = rfSetDevice.block(1, rfIndex, rfIndex);
+                Aurora::Matrix rfMf;
+                if(aPreComputes.matchedFilter.getDimSize(1) == 1) //Use CE
+                {
+                    rfMf = matchFilt(rf, Aurora::real(ifft(aPreComputes.matchedFilter.toDeviceMatrix())));
+                }
+                else // Use CEM
+                {
+                    rfMf = matchFilt(rf, Aurora::real(ifft(aPreComputes.matchedFilter.block(1, rfIndex, rfIndex).toDeviceMatrix())));
+                }
+                Aurora::Matrix rfWin = addWin(rfMf, transParams::offsetElectronic, transParams::minSpeedOfSound, transParams::maxSpeedOfSound, distance[cnt], transParams::aScanReconstructionFrequency);
+                Aurora::Matrix rfEnv = Aurora::abs(Aurora::hilbert(rfWin));
+                Aurora::Matrix rfNorm = rfEnv / Aurora::max(rfEnv);
+                tof[cnt] = findFirstvalueGreaterThanGivenValue(rfNorm, thre) + 1;
+            }
+        }
+    }
+    tof = offset + tof/transParams::aScanReconstructionFrequency;
+    Aurora::Matrix sosRough = distance / tof;
+    Aurora::Matrix mu = Aurora::mean(sosRough);
+    Aurora::Matrix sigma = Aurora::std(sosRough);
+    Aurora::Matrix idx_ = Aurora::abs(sosRough - mu) <= sigma;
+    Aurora::Matrix sosClean =  removeDataFromArrays(sosRough, idx_);
+    Aurora::Matrix sos = mean(sosClean);
+    RECON_INFO("SOS Value");
+    RECON_INFO(std::to_string(sos[0]));
+    return sos;
+}

src/common/SoundSpeedWaterCalculation/estimateSOSWater.h

@@ -8,7 +8,7 @@ class Parser;
 
 namespace Recon
 {
-    Aurora::Matrix estimateSOSWater(Parser* aParser, const CEInfo& aCE);
+    Aurora::Matrix estimateSOSWater(Parser* aParser, const CEInfo& aCE, const PreComputes& aPreComputes);
 }
 
 #endif

src/common/SoundSpeedWaterCalculation/getSoundSpeedWaterValue.cpp

@@ -0,0 +1,11 @@
+#include "getSoundSpeedWaterValue.h"
+#include "estimateSOSWater.h"
+
+using namespace Recon;
+
+SOSInfo Recon::getSoundSpeedWaterValue(Parser* aParser, const CEInfo& aCE, const PreComputes& aPreComputes)
+{
+    SOSInfo result;
+    result.expectedSOSWater = Recon::estimateSOSWater(aParser, aCE, aPreComputes);
+    return result;
+}

src/common/SoundSpeedWaterCalculation/getSoundSpeedWaterValue.h

@@ -0,0 +1,19 @@
+#ifndef GET_SOUND_SPEED_WATER_VALUE_H
+#define GET_SOUND_SPEED_WATER_VALUE_H
+
+#include "Matrix.h"
+#include "common/getMeasurementMetaData.h"
+
+class Parser;
+
+struct SOSInfo
+{
+    Aurora::Matrix expectedSOSWater;
+};
+
+namespace Recon
+{
+    SOSInfo getSoundSpeedWaterValue(Parser* aParser, const CEInfo& aCE, const PreComputes& aPreComputes);
+}
+
+#endif

src/common/SoundSpeedWaterCalculation/temperatureToSoundSpeed.cpp

@@ -0,0 +1,29 @@
+#include "temperatureToSoundSpeed.h"
+#include "estimateSOSWater.h"
+
+#include "Function1D.h"
+
+using namespace Recon;
+
+Aurora::Matrix Recon::temperatureToSoundSpeed(const Aurora::Matrix& aTemperature, const std::string& aMethod)
+{
+    Aurora::Matrix result;
+    if (aMethod == "marczak")
+    {
+        float* kData = new float[6] {2.78786e-9, -1.398845e-6, 3.287156e-4, -5.799136e-2, 5.038813, 1.402385e3};
+        Aurora::Matrix k = Aurora::Matrix::fromRawData(kData, 6);
+        result = polyval(k, aTemperature);
+    }
+    else if(aMethod == "mader")
+    {
+        float* kData = new float[6] {3.14643e-9, -1.478e-6, 0.000334199, -0.0580852, 5.03711, 1402.39};
+        Aurora::Matrix k = Aurora::Matrix::fromRawData(kData, 6);
+        result = polyval(k, aTemperature);
+    }
+    else if(aMethod == "jan")
+    {
+        float speedData = 1557 - 0.0245 * pow((74 - aTemperature.getData()[0]), 2);
+        result = Aurora::Matrix::fromRawData(new float[1] {speedData}, 1);
+    }
+    return result;
+}

src/common/SoundSpeedWaterCalculation/temperatureToSoundSpeed.h

@@ -0,0 +1,11 @@
+#ifndef TEMPERATURE_TO_SOUND_SPEED_H
+#define TEMPERATURE_TO_SOUND_SPEED_H
+
+#include "Matrix.h"
+
+namespace Recon
+{
+    Aurora::Matrix temperatureToSoundSpeed(const Aurora::Matrix& aTemperature, const std::string& aMethod);
+}
+
+#endif

src/common/getMeasurementMetaData.cpp

@@ -17,7 +17,7 @@
 #include "ceMatchedFilterHandling.h"
 #include "ShotList/ShotList.h"
 #include "config/config.h"
-#include "temperatureCalculation/estimateSOSWater.h"
+#include "common/SoundSpeedWaterCalculation/getSoundSpeedWaterValue.h"
 
 
 using namespace Recon;
@@ -127,36 +127,13 @@ TransFormInfo Recon::getTransformationMatrix(Parser* aParser, const Matrix& aMot
     return result;
 }
 
-Matrix Recon::temperatureToSoundSpeed(const Matrix& aTemperature, const std::string& aMethod)
+//该版本删除该函数
-{
+// TempInfo Recon::getTemperatureInfo(Parser* aParser, const CEInfo& aCEInfo, const PreComputes& aPreComputes)
-    Matrix result;
+// {
-    if (aMethod == "marczak")
+//     TempInfo result;
-    {
+//     result.expectedSOSWater = getSoundSpeedWaterValue(aParser,aCEInfo, aPreComputes);
-        float* kData = new float[6] {2.78786e-9, -1.398845e-6, 3.287156e-4, -5.799136e-2, 5.038813, 1.402385e3};
+//     return result;
-        Matrix k = Matrix::fromRawData(kData, 6);
+// }
-        result = polyval(k, aTemperature);
-    }
-    else if(aMethod == "mader")
-    {
-        float* kData = new float[6] {3.14643e-9, -1.478e-6, 0.000334199, -0.0580852, 5.03711, 1402.39};
-        Matrix k = Matrix::fromRawData(kData, 6);
-        result = polyval(k, aTemperature);
-    }
-    else if(aMethod == "jan")
-    {
-        float speedData = 1557 - 0.0245 * pow((74 - aTemperature.getData()[0]), 2);
-        result = Matrix::fromRawData(new float[1] {speedData}, 1);
-    }
-    return result;
-}
-
-//已验证，完全正确
-TempInfo Recon::getTemperatureInfo(Parser* aParser, const CEInfo& aCEInfo)
-{
-    TempInfo result;
-    result.expectedSOSWater = estimateSOSWater(aParser,aCEInfo);
-    return result;
-}
 
 //已验证 完全正确
 CEInfo Recon::getCEInfo(Parser* aParser, const MeasurementInfo aInfo)

src/common/getMeasurementMetaData.h

@@ -71,9 +71,8 @@ namespace Recon
     Aurora::Matrix getAvailableMotorPositions(Parser* aParser);
     MeasurementInfo loadMeasurementInfos(Parser* aParser);
     TransFormInfo getTransformationMatrix(Parser* aParser, const Aurora::Matrix& aMotorPosList);
-    TempInfo getTemperatureInfo(Parser* aParser, const CEInfo& aCEInfo);
+    //TempInfo getTemperatureInfo(Parser* aParser, const CEInfo& aCEInfo, const PreComputes& aPreComputes);
     CEInfo getCEInfo(Parser* aParser, const MeasurementInfo aInfo);
-    Aurora::Matrix temperatureToSoundSpeed(const Aurora::Matrix& aTemperature, const std::string& aMethod);
 
 }
 #endif // !GET_MEASUREMENT_METADATA_H

src/common/temperatureCalculation/estimateSOSWater.cpp

@@ -1,102 +0,0 @@
-#include "estimateSOSWater.h"
-#include "Function1D.cuh"
-#include "Function1D.h"
-#include "Matrix.h"
-#include "CudaMatrix.h"
-#include "common/dataBlockCreation/getAscanBlock.h"
-#include "transmissionReconstruction/detection/getTransmissionData.cuh"
-#include "transmissionReconstruction/detection/detection.h"
-#include "config/config.h"
-#include "config/geometryConfig.h"
-
-#include "Function2D.cuh"
-#include "Function2D.h"
-#include "common/dataBlockCreation/getAscanBlock.h"
-
-#include <cstring>
-
-using namespace Recon;
-
-namespace
-{
-    const Aurora::Matrix EMIT_TAS = Aurora::Matrix::fromRawData(new float[17]{94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110}, 17);
-    const Aurora::Matrix RECEIVE_TAS = Aurora::Matrix::fromRawData(new float[17]{103, 104, 105, 106, 107, 108, 109, 110, 94, 95, 96, 97, 98, 99, 100, 101, 102}, 17);
-    const Aurora::Matrix ALL_RECEIVER_TAS_LIST = Aurora::Matrix::fromRawData(new float[128] {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128},128); //1~128
-    const Aurora::Matrix ALL_RECEIVER_ELEMENT_LIST = Aurora::Matrix::fromRawData(new float[18] {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18},18);
-    const float SOS_INITIAL = 1509.3;
-    const float THRE = 0.5;
-    inline int findIndexFromEmitterAndReceiver(const Aurora::Matrix& aEmitter, float aEmitterValue, const Aurora::Matrix& aReceiver, float aReceiverValue)
-    {
-        for (int i = 0; i < aEmitter.getDataSize(); ++i)
-        {
-            if (aEmitter[i] == aEmitterValue && aReceiver[i] == aReceiverValue)
-            {
-                return i;
-            }
-        }
-        return -1;
-    }
-
-    Aurora::Matrix matchFilt(const Aurora::CudaMatrix& aSig, const Aurora::CudaMatrix& aMod)
-    {
-        Aurora::CudaMatrix fftMod = fft(aMod);
-        Aurora::CudaMatrix fftSig = fft(aSig);
-        Aurora::CudaMatrix mfFftSig = getTransmissionDataSubFunction(fftSig, fftMod);
-        Aurora::Matrix mfSig = Aurora::real(ifft(mfFftSig)).toHostMatrix();
-        std::memmove(mfSig.getData(), mfSig.getData() + 1, (mfSig.getDataSize() - 1) * sizeof(float));
-        return mfSig;
-    }
-
-    inline int findFirstvalueGreaterThanGivenValue(const Aurora::Matrix& aMatrix, float aValue)
-    {
-        for(int i=0; i<aMatrix.getDataSize(); ++i)
-        {
-            if(aMatrix[i] > aValue)
-            {
-                return i;
-            }
-        }
-        return -1;
-    }
-}
-
-Aurora::Matrix Recon::estimateSOSWater(Parser *aParser, const CEInfo &aCE)
-{
-    int offset = transParams::priorWvOffset;
-    int rn = 18;
-    int emitterSize = EMIT_TAS.getDataSize();
-    Aurora::Matrix distance = Aurora::Matrix::fromRawData(new float[emitterSize]{0}, emitterSize);
-    Aurora::Matrix mp = Aurora::Matrix::fromRawData(new float[1]{1}, 1);
-    Aurora::Matrix sn = Aurora::Matrix::fromRawData(new float[1]{18}, 1);
-    Aurora::Matrix tof = Aurora::Matrix::fromRawData(new float[emitterSize], 1, emitterSize);
-    float snValue = 18;
-    #pragma omp parallel for
-    for(int i=0; i<emitterSize; ++i)
-    {
-        Aurora::Matrix sl = Aurora::Matrix::fromRawData(new float[1]{EMIT_TAS[i]}, 1);
-        float slValue = EMIT_TAS[i];
-        float rl = RECEIVE_TAS[i];
-        Aurora::Matrix positionOfEmitter = emitterPositions[slValue-1](snValue-1).toMatrix();
-        Aurora::Matrix positionOfReceiver = emitterPositions[rl-1](rn-1).toMatrix();
-        distance[i] = Aurora::norm(positionOfEmitter - positionOfReceiver, Aurora::Norm2);
-        AscanBlock ascan = getAscanBlock(aParser, mp, sl, sn, ALL_RECEIVER_TAS_LIST, ALL_RECEIVER_ELEMENT_LIST);
-        size_t dataSize = ascan.ascanBlock.getDataSize();
-        short* data = new short[dataSize];
-        std::copy(ascan.ascanBlock.getData(), ascan.ascanBlock.getData() + dataSize, data);
-        Aurora::Matrix rf = Aurora::convertfp16tofloat(data, ascan.ascanBlock.getDimSize(0), ascan.ascanBlock.getDimSize(1));
-        delete[] data;
-        Aurora::CudaMatrix rfDevice = rf.toDeviceMatrix() / ascan.gainBlock.toDeviceMatrix();
-        int rfIndex = findIndexFromEmitterAndReceiver(aCE.tasIndices, rl, aCE.receiverIndices, rn);
-        rfDevice = rfDevice.block(1, rfIndex, rfIndex);
-        Aurora::Matrix rfMf = matchFilt(rfDevice, aCE.ceRef.toDeviceMatrix());
-        Aurora::Matrix rfWin = applyTimeWindowing(rfMf, transParams::aScanReconstructionFrequency, distance(i).toMatrix(), SOS_INITIAL, 
-                                                  transParams::offsetElectronic * transParams::aScanReconstructionFrequency, transParams::detectionWindowSOS,
-                                                  transParams::minSpeedOfSound, transParams::maxSpeedOfSound, transParams::gaussWindow).AscanBlockProcessed;
-        Aurora::Matrix rfEnv = abs(Aurora::hilbert(rfWin));
-        Aurora::Matrix rfNorm = rfEnv / max(rfEnv);
-        tof[i] = findFirstvalueGreaterThanGivenValue(rfNorm, THRE) + 1;
-    }
-    tof = (tof + offset) / transParams::aScanReconstructionFrequency;
-    auto sos = mean(distance / tof);
-    return sos;
-}

src/reflectionReconstruction/preprocessData/preprocessTransmissionReconstructionForReflection.cpp

@@ -18,7 +18,7 @@ using namespace Recon;
 
 PreprocessReflectionData Recon::preprocessTransmissionReconstructionForReflection(const Aurora::Matrix& aSOSMap, const Aurora::Matrix& aATTMap,
                                                          const Aurora::Matrix& aDdims, const GeometryInfo& aGeometryInfo,
-                                                         const TempInfo& aTemp)
+                                                         const SOSInfo& aSOSInfo)
 {
     PreprocessReflectionData result;
     int soundSpeedCorrection = reflectParams::soundSpeedCorrection;
@@ -57,7 +57,7 @@ PreprocessReflectionData Recon::preprocessTransmissionReconstructionForReflectio
         {
             if(std::abs(sppedMap3d[i]) < 1)
             {
-                sppedMap3d[i] = aTemp.expectedSOSWater[0];
+                sppedMap3d[i] = aSOSInfo.expectedSOSWater[0];
             }
         }
         result.transRecos.speedMap3D = sppedMap3d;
@@ -116,7 +116,7 @@ PreprocessReflectionData Recon::preprocessTransmissionReconstructionForReflectio
         //     transRecos.SpeedMap3D = reshape(temperatureToSoundSpeed(polyvaln(temp.TemperatureModel4D.p,[XNEW(:),YNEW(:), ZNEW(:), zeros(size(XNEW(:)))])),size(XNEW));
         // catch
 
-        result.transRecos.speedMap3D = zeros(x.getDimSize(1), y.getDimSize(1), z.getDimSize(1)) + aTemp.expectedSOSWater;
+        result.transRecos.speedMap3D = zeros(x.getDimSize(1), y.getDimSize(1), z.getDimSize(1)) + aSOSInfo.expectedSOSWater;
     
         
         result.transRecos.beginTransMap = begin_TransMap;

src/reflectionReconstruction/preprocessData/preprocessTransmissionReconstructionForReflection.h

@@ -4,6 +4,7 @@
 #include "Matrix.h"
 #include "common/getGeometryInfo.h"
 #include "common/getMeasurementMetaData.h"
+#include "common/SoundSpeedWaterCalculation/getSoundSpeedWaterValue.h"
 
 namespace Recon
 {
@@ -24,7 +25,7 @@ namespace Recon
     };
     PreprocessReflectionData preprocessTransmissionReconstructionForReflection(const Aurora::Matrix& aSOSMap, const Aurora::Matrix& aATTMap,
                                                              const Aurora::Matrix& aDdims, const GeometryInfo& aGeometryInfo,
-                                                             const TempInfo& aTemp);
+                                                             const SOSInfo& aSOSInfo);
 }
 
 #endif

src/startReconstructions.cpp

@@ -20,9 +20,10 @@
 #include "common/estimatePulseLength.h"
 #include "common/ceMatchedFilterHandling.h"
 #include "common/createPositioningImage.h"
+#include "common/estimateOffset.h"
+#include "common/SoundSpeedWaterCalculation/getSoundSpeedWaterValue.h"
 #include "transmissionReconstruction/startTransmissionReconstruction.h"
 #include "transmissionReconstruction/saveTransmissionImagesInReflCoords.h"
-#include "reflectionReconstruction/preprocessData/estimateOffset.h"
 #include "reflectionReconstruction/preprocessData/precalcImageParameters.h"
 #include "reflectionReconstruction/preprocessData/preprocessTransmissionReconstructionForReflection.h"
 #include "reflectionReconstruction/startReflectionReconstruction.h"
@@ -86,21 +87,38 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
     float maxNumTAS = max(auroraUnion(slList, rlList)).getData()[0];
     MeasurementInfo expInfo = loadMeasurementInfos(&dataParser);
     CEInfo ce = getCEInfo(&dataParser, expInfo);
-    TempInfo temp = getTemperatureInfo(&dataParser, ce);
+
+    PreComputes preComputes;
+    if(ce.measuredCEUsed)
+    {
+        preComputes.matchedFilter = createMatchedFilter(ce.ce, ce.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
+    }
+    else
+    {
+        preComputes.matchedFilter = createMatchedFilter(ce.ceRef, ce.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
+    }
+    preComputes.timeInterval = (float)1 / reflectParams::aScanReconstructionFrequency;
+    preComputes.measuredCEUsed = ce.measuredCEUsed;
+    preComputes.measuredCE_TASIndices = ce.tasIndices;
+    preComputes.measuredCE_receiverIndices = ce.receiverIndices;
+    preComputes.offset = estimateOffset(expInfo, ce);
+
+    SOSInfo sosValue = getSoundSpeedWaterValue(&dataParser, ce, preComputes);
+    SOSInfo sosValueRef;
+
     TransFormInfo transformationInfo = getTransformationMatrix(&dataParser, motorPosTotal);
     Matrix transformationMatrices = transformationInfo.rotationMatrix;
     Matrix motorPosAvailable = transformationInfo.motorPos;
 
     Matrix motorPosAvailableRef;
     MeasurementInfo expInfoRef;
-    TempInfo tempRef;
     CEInfo ceRef;
     Matrix transformationMatricesRef;
     if(transParams::runTransmissionReco)
     {
         expInfoRef = loadMeasurementInfos(&refParser);
         ceRef = getCEInfo(&refParser, expInfoRef);
-        tempRef = getTemperatureInfo(&refParser, ceRef);
+        sosValueRef = getSoundSpeedWaterValue(&refParser, ceRef, preComputes);
         transformationInfo = getTransformationMatrix(&refParser, motorPosTotal);
         transformationMatricesRef = transformationInfo.rotationMatrix;
         motorPosAvailableRef = transformationInfo.motorPos;
@@ -142,18 +160,6 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
     }
 
     GeometryInfo geom = getGeometryInfo(motorPosAvailable, transformationMatrices, rlList, rnList, slList, snList);
-    PreComputes preComputes;
-    if((reflectParams::matchedFilterCeAScan && reflectParams::runReflectionReco) || transParams::runTransmissionReco)
-    {
-        if(ce.measuredCEUsed)
-        {
-            preComputes.matchedFilter = createMatchedFilter(ce.ce, ce.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
-        }
-        else
-        {
-            preComputes.matchedFilter = createMatchedFilter(ce.ceRef, ce.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
-        }
-    }
     if(expInfo.sampleRate != reflectParams::aScanReconstructionFrequency)
     {
         reflectParams::expectedAScanDataLength = ceil(expInfo.numberSamples * ((float)reflectParams::aScanReconstructionFrequency / expInfo.sampleRate));
@@ -172,15 +178,8 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
         {
             preComputes.matchedFilterRef = createMatchedFilter(ceRef.ceRef, ceRef.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
         }
-        preComputes.timeInterval = (float)1 / reflectParams::aScanReconstructionFrequency;
-        preComputes.measuredCEUsed = ce.measuredCEUsed;
-        preComputes.measuredCE_TASIndices = ce.tasIndices;
-        preComputes.measuredCE_receiverIndices = ce.receiverIndices;
-        preComputes.offset = estimateOffset(expInfo, ce);
-
         expInfo.matchedFilter = preComputes.matchedFilter;
         expInfoRef.matchedFilter = preComputes.matchedFilterRef;
-
         CeEstimatePulseLength ceEstimatePulseLength;
         CeEstimatePulseLength ceEstimatePulseLengthRef;
         if(ce.measuredCEUsed)
@@ -217,7 +216,7 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
 
         GeometryInfo geomRef = getGeometryInfo(motorPosAvailableRef, transformationMatricesRef, rlList, rnList, slList, snList);
         RECON_INFO("Start transmissionRecostruction.");
-        transmissionResult = startTransmissionReconstruction(mp_inter, mpRef_inter, slList_inter, snList_inter, rlList_inter, rnList_inter, temp, tempRef, geom, geomRef, expInfo, expInfoRef, preComputes, &dataParser, &refParser);
+        transmissionResult = startTransmissionReconstruction(mp_inter, mpRef_inter, slList_inter, snList_inter, rlList_inter, rnList_inter, sosValue, sosValueRef, geom, geomRef, expInfo, expInfoRef, preComputes, &dataParser, &refParser);
         attAvailable = true;
         sosAvailable = true;
         Matrix minPos =  Matrix::fromRawData(new float[3],3,1,1);
@@ -238,17 +237,12 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
         Matrix recoATT = transmissionResult.recoATT;      
         //检测可使用内存是否足够,输出警报用,todo
         //checkEnvAndMemory(reflectParams.imageInfos.IMAGE_XYZ);
-        auto preProcessData = preprocessTransmissionReconstructionForReflection(recoSOS, recoATT, transmissionResult.ddmis, geom, temp);
+        auto preProcessData = preprocessTransmissionReconstructionForReflection(recoSOS, recoATT, transmissionResult.ddmis, geom, sosValue);
         Matrix mp_inter = intersect(motorPosAvailable, reflectParams::motorPos);
         Matrix slList_inter = intersect(slList, reflectParams::senderTasList);
         Matrix snList_inter = intersect(snList, reflectParams::senderElementList);
         Matrix rlList_inter = intersect(rlList, reflectParams::receiverTasList);
         Matrix rnList_inter = intersect(rnList, reflectParams::receiverElementList);
-        preComputes.timeInterval = (float)1 / reflectParams::aScanReconstructionFrequency;
-        preComputes.measuredCEUsed = ce.measuredCEUsed;
-        preComputes.measuredCE_TASIndices = ce.tasIndices;
-        preComputes.measuredCE_receiverIndices = ce.receiverIndices;
-        preComputes.offset = estimateOffset(expInfo, ce);
 
         reflectParams::gpuSelectionList = reconParams::gpuSelectionList;
         RECON_INFO("Start reflectionRecostruction.");

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

@@ -78,7 +78,7 @@ namespace
     }
 
     BlockOfTransmissionData getBlockOfTransmissionData(const Matrix& aMp, const Matrix& aMpRef, const Matrix& aSl, const Matrix& aSn, const Matrix& aRlList, const Matrix& aRnList,
-                                    const TempInfo& aTemp, const TempInfo& aTempRef, GeometryInfo aGeom, GeometryInfo& aGeomRef,
+                                    const SOSInfo& aSOS, const SOSInfo& aSOSRef, GeometryInfo aGeom, GeometryInfo& aGeomRef,
                                     const Matrix& aSnrRmsNoise, const Matrix& aSnrRmsNoiseRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
                                     const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
     {
@@ -162,13 +162,13 @@ namespace
 }
 
 void  getBlockOfTransmissionDataInThread(size_t aIndex, const Matrix& aMp, const Matrix& aMpRef, const Matrix& aSl, const Matrix& aSn, const Matrix& aRlList, const Matrix& aRnList,
-                                         const TempInfo& aTemp, const TempInfo& aTempRef, GeometryInfo aGeom, GeometryInfo aGeomRef,
+                                         const SOSInfo& aSOS, const SOSInfo& aSOSRef, GeometryInfo aGeom, GeometryInfo aGeomRef,
                                          const Matrix& aSnrRmsNoise, const Matrix& aSnrRmsNoiseRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
                                          const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef, unsigned int aGPUId)
 {
     cudaSetDevice(aGPUId); 
-    auto buffer = getBlockOfTransmissionData(aMp, aMpRef, aSl, aSn, aRlList, aRnList, aTemp, 
+    auto buffer = getBlockOfTransmissionData(aMp, aMpRef, aSl, aSn, aRlList, aRnList, aSOS, 
-                                             aTempRef, aGeom, aGeomRef, aSnrRmsNoise, aSnrRmsNoiseRef,
+                                             aSOSRef, aGeom, aGeomRef, aSnrRmsNoise, aSnrRmsNoiseRef,
                                              aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef);
     std::unique_lock<std::mutex> lock(CREATE_BUFFER_MUTEX); 
     BLOCK_OF_TRANSIMISSIONDARA_BUFFER[std::to_string(aIndex)] = buffer;
@@ -178,7 +178,7 @@ void  getBlockOfTransmissionDataInThread(size_t aIndex, const Matrix& aMp, const
 }
 
 void  createThreadForGetBlockOfTransmissionData(const Matrix& aMotorPos, const Matrix& aMotoPosRef, const Matrix& aSlList, const Matrix& aSnList, const Matrix& aRlList, const Matrix& aRnList,
-                                         const TempInfo& aTemp, const TempInfo& aTempRef, GeometryInfo aGeom, GeometryInfo aGeomRef,
+                                         const SOSInfo& aSOS, const SOSInfo& aSOSRef, GeometryInfo aGeom, GeometryInfo aGeomRef,
                                          const Matrix& aSnrRmsNoise, const Matrix& aSnrRmsNoiseRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
                                          const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
 {
@@ -202,7 +202,7 @@ void  createThreadForGetBlockOfTransmissionData(const Matrix& aMotorPos, const M
                 CREATE_BUFFER_CONDITION.wait(lock, []{return BUFFER_COUNT<BUFFER_SIZE;});
                 ++BUFFER_COUNT;
                 lock.unlock();
-                speedUpThread[index] =  std::thread(getBlockOfTransmissionDataInThread,index,mp,mpRef,sl,sn,aRlList,aRnList,aTemp,aTempRef,aGeom,aGeomRef,aSnrRmsNoise,aSnrRmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef, index % BUFFER_SIZE);
+                speedUpThread[index] =  std::thread(getBlockOfTransmissionDataInThread,index,mp,mpRef,sl,sn,aRlList,aRnList,aSOS,aSOSRef,aGeom,aGeomRef,aSnrRmsNoise,aSnrRmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef, index % BUFFER_SIZE);
             }
         }
     }
@@ -216,7 +216,7 @@ void  createThreadForGetBlockOfTransmissionData(const Matrix& aMotorPos, const M
 
 TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, const Aurora::Matrix& aMotoPosRef, const Aurora::Matrix& aSlList,
                                 const Aurora::Matrix& aSnList, const Aurora::Matrix& aRlList, const Aurora::Matrix& aRnList,
-                                const TempInfo& aTemp, const TempInfo& aTempRef, GeometryInfo& aGeom,
+                                const SOSInfo& aSOS, const SOSInfo& aSOSRef, GeometryInfo& aGeom,
                                 GeometryInfo& aGeomRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
                                 const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
 {
@@ -259,7 +259,7 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
         rnBlockTotal = zeros(1,numScans,1);
     }
 
-    std::thread speedUpThread = std::thread(createThreadForGetBlockOfTransmissionData,aMotorPos,aMotoPosRef,aSlList,aSnList,aRlList,aRnList,aTemp,aTempRef,aGeom,aGeomRef,rmsNoise,rmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef);
+    std::thread speedUpThread = std::thread(createThreadForGetBlockOfTransmissionData,aMotorPos,aMotoPosRef,aSlList,aSnList,aRlList,aRnList,aSOS,aSOSRef,aGeom,aGeomRef,rmsNoise,rmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef);
     int numData = 0;
     int numPossibleScans = 0;
     sched_param sch;
@@ -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(), 
-                                                             aTemp.expectedSOSWater[0], aTempRef.expectedSOSWater[0]);
+                                                             aSOS.expectedSOSWater[0], aSOSRef.expectedSOSWater[0], aPreComputes.offset);
                 blockData.attData = detect.att;
                 blockData.tofData = detect.tof;
                 BlockOfTransmissionData transmissionBlock=blockData;

src/transmissionReconstruction/detection/getTransmissionData.h

@@ -4,6 +4,7 @@
 #include "Matrix.h"
 #include "common/getMeasurementMetaData.h"
 #include "common/getGeometryInfo.h"
+#include "common/SoundSpeedWaterCalculation/getSoundSpeedWaterValue.h"
 
 class Parser;
 namespace Recon
@@ -42,7 +43,7 @@ namespace Recon
 
     TransmissionData getTransmissionData(const Aurora::Matrix& aMotorPos, const Aurora::Matrix& aMotoPosRef, const Aurora::Matrix& aSlList,
                              const Aurora::Matrix& aSnList, const Aurora::Matrix& aRlList, const Aurora::Matrix& aRnList,
-                             const TempInfo& aTemp, const TempInfo& aTempRef, GeometryInfo& aGeom,
+                             const SOSInfo& aSOS, const SOSInfo& aSOSRef, GeometryInfo& aGeom,
                              GeometryInfo& aGeomRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
                              const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef);
 }

src/transmissionReconstruction/startTransmissionReconstruction.cpp

@@ -19,16 +19,16 @@ using namespace Recon;
 
 TransmissionReconstructionResult Recon::startTransmissionReconstruction(const Aurora::Matrix& aMotorPos, const Aurora::Matrix& aMotoPosRef, const Aurora::Matrix& aSlList,
                                      const Aurora::Matrix& aSnList, const Aurora::Matrix& aRlList, const Aurora::Matrix& aRnList,
-                                     const TempInfo& aTemp, const TempInfo& aTempRef, Recon::GeometryInfo& aGeom,
+                                     const SOSInfo& aSOS, const SOSInfo& aSOSRef, Recon::GeometryInfo& aGeom,
                                      GeometryInfo& aGeomRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
                                      const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
 {
     RECON_INFO("Start getTransmissionData.");
-    auto transmissionData = getTransmissionData(aMotorPos, aMotoPosRef, aSlList, aSnList, aRlList, aRnList, aTemp, aTempRef, 
+    auto transmissionData = getTransmissionData(aMotorPos, aMotoPosRef, aSlList, aSnList, aRlList, aRnList, aSOS, aSOSRef, 
                                                 aGeom, aGeomRef, aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef);    
     Matrix dists = Recon::distanceBetweenTwoPoints(transmissionData.senderList, transmissionData.receiverList);
 
-    Matrix valid = Recon::checkTofDetections(transmissionData.tofDataTotal, dists, aTemp.expectedSOSWater[0], 
+    Matrix valid = Recon::checkTofDetections(transmissionData.tofDataTotal, dists, aSOS.expectedSOSWater[0], 
                                              Recon::transParams::minSpeedOfSound,Recon::transParams::maxSpeedOfSound).valid;    
     if(transParams::qualityCheck)
     {
@@ -40,7 +40,7 @@ TransmissionReconstructionResult Recon::startTransmissionReconstruction(const Au
     Matrix senderList = removeDataFromArrays(positionValues.senderCoordList, valid);
     Matrix reveiverList = removeDataFromArrays(positionValues.receiverCoordList, valid);
     RECON_INFO("Start reconstructArt.");
-    auto transmissionReon = reconstructArt(tofData, attData, positionValues.dims, senderList, reveiverList, positionValues.res, aTemp.expectedSOSWater[0]);
+    auto transmissionReon = reconstructArt(tofData, attData, positionValues.dims, senderList, reveiverList, positionValues.res, aSOS.expectedSOSWater[0]);
 
     TransmissionReconstructionResult result;
     result.recoATT = transmissionReon.outATT;

src/transmissionReconstruction/startTransmissionReconstruction.h

@@ -5,6 +5,7 @@
 #include "common/getMeasurementMetaData.h"
 #include "common/getGeometryInfo.h"
 #include "transmissionReconstruction/detection/getTransmissionData.h"
+#include "common/SoundSpeedWaterCalculation/getSoundSpeedWaterValue.h"
 
 class Parser;
 
@@ -19,7 +20,7 @@ namespace Recon
 
     TransmissionReconstructionResult startTransmissionReconstruction(const Aurora::Matrix& aMotorPos, const Aurora::Matrix& aMotoPosRef, const Aurora::Matrix& aSlList,
                                          const Aurora::Matrix& aSnList, const Aurora::Matrix& aRlList, const Aurora::Matrix& aRnList, 
-                                         const TempInfo& aTemp, const TempInfo& aTempRef, Recon::GeometryInfo& aGeom,
+                                         const SOSInfo& aTemp, const SOSInfo& aTempRef, Recon::GeometryInfo& aGeom,
                                          GeometryInfo& aGeomRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
                                          const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef); 
 }