diff --git a/src/common/getMeasurementMetaData.cpp b/src/common/getMeasurementMetaData.cpp index df9ab3f..c21cf59 100644 --- a/src/common/getMeasurementMetaData.cpp +++ b/src/common/getMeasurementMetaData.cpp @@ -6,6 +6,7 @@ #include "Function.h" #include "Function2D.h" +#include "Matrix.h" #include "Parser.h" #include "Data/MetaData.h" #include "Data/MovementData.h" @@ -16,6 +17,7 @@ #include "ceMatchedFilterHandling.h" #include "ShotList/ShotList.h" #include "config/config.h" +#include "temperatureCalculation/estimateSOSWater.h" using namespace Recon; @@ -148,110 +150,10 @@ Matrix Recon::temperatureToSoundSpeed(const Matrix& aTemperature, const std::str } //已验证,完全正确 -TempInfo Recon::getTemperatureInfo(Parser* aParser, float aNumTas) +TempInfo Recon::getTemperatureInfo(Parser* aParser, const CEInfo& aCEInfo) { TempInfo result; - //jumoTemp - JumoTemperaturePointer jumoTemp1 = aParser->getTemperatureData().getJumoTemperature1(); - JumoTemperaturePointer jumoTemp2 = aParser->getTemperatureData().getJumoTemperature1(); - JumoTemperaturePointer jumoTemp3 = aParser->getTemperatureData().getJumoTemperature1(); - JumoTemperaturePointer jumoTemp4 = aParser->getTemperatureData().getJumoTemperature1(); - int jumoTempNum = jumoTemp1.getLength(); - if(jumoTemp2.getLength() < jumoTempNum) - { - jumoTempNum = jumoTemp2.getLength(); - } - - if(jumoTemp3.getLength() < jumoTempNum) - { - jumoTempNum = jumoTemp3.getLength(); - } - - if(jumoTemp4.getLength() < jumoTempNum) - { - jumoTempNum = jumoTemp4.getLength(); - } - - float* jumoTempData = Aurora::malloc(jumoTempNum * 4); - for(int i=0; igetTemperatureData().getTasTemperature()[0].getLength(); - float* fromTasTempData = aParser->getTemperatureData().getTasTemperature()[0].get(); - float* tasTempData = new float[tasTempLength]; - std::copy(fromTasTempData, fromTasTempData+ tasTempLength, tasTempData); - Matrix tasTemp = Matrix::fromRawData(tasTempData, 2, tasTempLength / 2); - - if (!reconParams::correctTASTemp || result.jumoTemp.isNan()) - { - result.tasTemperature = tasTemp; - } - else - { - //todo JumoTemp都没有值 - // temp.TASTemperature(1, :, :) = correctTASTemperatures(temp2, temp.jumoTemp); - - // catch ME - - // % if tas temperatures not there - // % take temperatures from calibrated sensors - // % and use them for all TAS the same - - // writeReconstructionLog(sprintf('%s: %s TAS temperatures approximated from calibrated reference sensors.', ME.identifier, ME.message), 3) - - // if(~jumoAllNaN) - // temp.TASTemperature = zeros(2, numTAS, size(temp.jumoTemp, 2)); - // temp.TASTemperature(1, :, :) = repmat(mean(temp.jumoTemp), size(temp.TASTemperature, 2), 1); - // temp.TASTemperature(2, :, :) = repmat([1:size(temp.TASTemperature, 2)]', 1, size(temp.TASTemperature, 3)); - // else % worst case: no temperature data at all! cannot continue! - // error([mfilename ':noTemperatureAvailable'],'No TAS and no Jumo temperature available! Cannot reconstruct!'); - // end - - // end - } - } - - if (result.jumoTemp.isNan()) - { - result.expectedTemp = mean(result.tasTemperature(0,$,$).toMatrix(),FunctionDirection::Column,false); - result.jumoTemp = result.expectedTemp; - } - else - { - result.expectedTemp = mean(result.jumoTemp,FunctionDirection::All,false); - } - result.expectedSOSWater = temperatureToSoundSpeed(result.expectedTemp , "marczak"); - + result.expectedSOSWater = estimateSOSWater(aParser,aCEInfo); return result; } diff --git a/src/common/getMeasurementMetaData.h b/src/common/getMeasurementMetaData.h index ff1f17e..be56d26 100644 --- a/src/common/getMeasurementMetaData.h +++ b/src/common/getMeasurementMetaData.h @@ -44,9 +44,6 @@ namespace Recon struct TempInfo { - Aurora::Matrix jumoTemp; - Aurora::Matrix tasTemperature; - Aurora::Matrix expectedTemp; Aurora::Matrix expectedSOSWater; }; @@ -72,7 +69,7 @@ namespace Recon Aurora::Matrix getAvailableMotorPositions(Parser* aParser); MeasurementInfo loadMeasurementInfos(Parser* aParser); TransFormInfo getTransformationMatrix(Parser* aParser, const Aurora::Matrix& aMotorPosList); - TempInfo getTemperatureInfo(Parser* aParser, float aNumTas); + TempInfo getTemperatureInfo(Parser* aParser, const CEInfo& aCEInfo); CEInfo getCEInfo(Parser* aParser, const MeasurementInfo aInfo); Aurora::Matrix temperatureToSoundSpeed(const Aurora::Matrix& aTemperature, const std::string& aMethod); diff --git a/src/common/temperatureCalculation/estimateSOSWater.cpp b/src/common/temperatureCalculation/estimateSOSWater.cpp new file mode 100644 index 0000000..ec5d801 --- /dev/null +++ b/src/common/temperatureCalculation/estimateSOSWater.cpp @@ -0,0 +1,102 @@ +#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 + +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 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(); - } + } + if(detection.contains("priorWvOffset")) + { + transParams::priorWvOffset = detection.at("priorWvOffset").get(); + } } if(params.contains("rayTracing")) @@ -748,6 +752,7 @@ namespace Recon transParams::detectionWindowATT = 50; transParams::pulseLengthSamples = 0; transParams::pulseLengthRefSamples = 0; + transParams::priorWvOffset = -3; //transParams.rayTracing transParams::bentReconstruction = false; transParams::bresenham = 1; diff --git a/src/config/config.h b/src/config/config.h index 10ab7d5..7dbba7c 100644 --- a/src/config/config.h +++ b/src/config/config.h @@ -146,6 +146,7 @@ namespace Recon EXTERN_C int detectionWindowATT; EXTERN_C float pulseLengthSamples; EXTERN_C float pulseLengthRefSamples; + EXTERN_C int priorWvOffset; //transParams.rayTracing EXTERN_C bool bentReconstruction; EXTERN_C int bresenham; diff --git a/src/startReconstructions.cpp b/src/startReconstructions.cpp index ced7469..db0b974 100644 --- a/src/startReconstructions.cpp +++ b/src/startReconstructions.cpp @@ -83,8 +83,8 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string //getMeasurementMetaData float maxNumTAS = max(auroraUnion(slList, rlList)).getData()[0]; MeasurementInfo expInfo = loadMeasurementInfos(&dataParser); - TempInfo temp = getTemperatureInfo(&dataParser, maxNumTAS); CEInfo ce = getCEInfo(&dataParser, expInfo); + TempInfo temp = getTemperatureInfo(&dataParser, ce); TransFormInfo transformationInfo = getTransformationMatrix(&dataParser, motorPosTotal); Matrix transformationMatrices = transformationInfo.rotationMatrix; Matrix motorPosAvailable = transformationInfo.motorPos; @@ -97,8 +97,8 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string if(transParams::runTransmissionReco) { expInfoRef = loadMeasurementInfos(&refParser); - tempRef = getTemperatureInfo(&refParser, maxNumTAS); ceRef = getCEInfo(&refParser, expInfoRef); + tempRef = getTemperatureInfo(&refParser, ceRef); transformationInfo = getTransformationMatrix(&refParser, motorPosTotal); transformationMatricesRef = transformationInfo.rotationMatrix; motorPosAvailableRef = transformationInfo.motorPos; diff --git a/src/transmissionReconstruction/dataFilter/dataFilter.cpp b/src/transmissionReconstruction/dataFilter/dataFilter.cpp index 9530916..be76a89 100644 --- a/src/transmissionReconstruction/dataFilter/dataFilter.cpp +++ b/src/transmissionReconstruction/dataFilter/dataFilter.cpp @@ -75,7 +75,7 @@ namespace Recon { } checkTofDetectionsResult checkTofDetections(const Aurora::Matrix &aVTofValues, const Aurora::Matrix &aVDists, - const Aurora::Matrix &aVSosRef, + float aVSosRef, float minSpeedOfSound, float maxSpeedOfSound) { @@ -122,11 +122,10 @@ namespace Recon { for (size_t i = 0; i < aVSNRList.getDataSize(); i++) { if (finite_SNRList.getData()[i] == 1.0){ - meanSNR+=aVSNRList[i]; std_SNRListData[j++] = aVSNRList[i]; } } - meanSNR = meanSNR/(float)j; + meanSNR = mean(aVSNRList)[0]; Aurora::Matrix std_SNRList = Aurora::Matrix::New(std_SNRListData,count,1,1); std_SNRList = Aurora::std(std_SNRList); float localSNRDifference = 2 * std_SNRList.getScalar(); diff --git a/src/transmissionReconstruction/dataFilter/dataFilter.h b/src/transmissionReconstruction/dataFilter/dataFilter.h index 7ee9b76..1f10568 100644 --- a/src/transmissionReconstruction/dataFilter/dataFilter.h +++ b/src/transmissionReconstruction/dataFilter/dataFilter.h @@ -27,7 +27,7 @@ struct checkTofDetectionsResult checkTofDetectionsResult checkTofDetections(const Aurora::Matrix &aVTofValues, const Aurora::Matrix &aVDists, - const Aurora::Matrix &aVSosRef, + float aVSosRef, float minSpeedOfSound, float maxSpeedOfSound); /** diff --git a/src/transmissionReconstruction/detection/detection.cpp b/src/transmissionReconstruction/detection/detection.cpp index 10170d1..f1cbc12 100644 --- a/src/transmissionReconstruction/detection/detection.cpp +++ b/src/transmissionReconstruction/detection/detection.cpp @@ -8,6 +8,7 @@ #include "Function2D.h" #include "Function3D.h" +#include "Matrix.h" #include "common/getMeasurementMetaData.h" #include "config/config.h" #include "transmissionReconstruction/detection/detection.cuh" @@ -81,19 +82,18 @@ namespace Recon { } TimeWindowResult applyTimeWindowing(const Aurora::Matrix &AscanBlock, float sampleRate, - const Aurora::Matrix &distBlock, const Aurora::Matrix &sosBlock, - float expectedSOSWater, float startOffset, float segmentLenOffset, + const Aurora::Matrix &distBlock, float sosWater, + float startOffset, float segmentLenOffset, float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow) { - auto sosOffset = calculateSOSOffset(sosBlock, expectedSOSWater, distBlock, sampleRate); - + Aurora::Matrix sosOffset = Aurora::Matrix::fromRawData(new float[1]{0}, 1); auto calcResult = calculateStarEndSearchPosition(distBlock, minSpeedOfSound, maxSpeedOfSound, sampleRate,AscanBlock.getDimSize(0), sosOffset, startOffset, segmentLenOffset); auto AscanBlockProcessed = zeros(AscanBlock.getDimSize(0),AscanBlock.getDimSize(1)); if(gaussWindow) { - auto expectedPosWater = (distBlock / expectedSOSWater) * sampleRate + startOffset; + auto expectedPosWater = (distBlock / sosWater) * sampleRate + startOffset; auto windowWidth = calcResult.endSearch-calcResult.startSearch; #pragma omp parallel for for (size_t i = 0; i < AscanBlock.getDimSize(1); i++) @@ -132,7 +132,7 @@ namespace Recon { Aurora::Matrix detectAttVectorized(const Aurora::Matrix &Ascan, const Aurora::Matrix &AscanRef, const Aurora::Matrix &distRef, - const Aurora::Matrix &sosWaterRef, + float sosWaterRef, const Aurora::Matrix &tof, int aScanReconstructionFrequency, float offsetElectronic, int detectionWindowATT) { @@ -173,8 +173,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 expectedSOSWater, + float aSOSWater, float aSOSWaterRef, int useTimeWindowing, int aScanReconstructionFrequency, float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow) @@ -188,12 +187,12 @@ namespace Recon { timeResult2.AscanBlockProcessed = AscanRefBlock; if (useTimeWindowing == 1) { timeResult1 = applyTimeWindowing( - AscanBlock, sampleRate, distBlock, sosWaterBlock, - expectedSOSWater, offsetElectronicSamples, detectionWindowSOS, + AscanBlock, sampleRate, distBlock, + aSOSWater, offsetElectronicSamples, detectionWindowSOS, minSpeedOfSound, maxSpeedOfSound, gaussWindow); timeResult2 = applyTimeWindowing( - AscanRefBlock, sampleRate, distBlockRef, sosWaterRefBlock, - expectedSOSWater, offsetElectronicSamples, detectionWindowSOS, + AscanRefBlock, sampleRate, distBlockRef, + aSOSWaterRef, offsetElectronicSamples, detectionWindowSOS, minSpeedOfSound, maxSpeedOfSound, gaussWindow); diffStartSearch = timeResult1.startSearch - timeResult2.startSearch; @@ -234,7 +233,7 @@ namespace Recon { if (useTimeWindowing) { shiftInSamples = shiftInSamples - diffStartSearch; } - auto tof = shiftInSamples / sampleRate + distBlock / sosWaterBlock; + auto tof = shiftInSamples / sampleRate + distBlock / aSOSWater; auto sosValue = distBlock / tof; DetectResult result; result.tof = tof; @@ -245,9 +244,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 expectedSOSWater, + int resampleFactor,int nthreads, float aSOSWater, float aSOSWaterRef, int useTimeWindowing, int aScanReconstructionFrequency,int detectionWindowATT, float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow) @@ -261,12 +258,12 @@ namespace Recon { timeResult2.AscanBlockProcessed = AscanRefBlock; if (useTimeWindowing == 1) { timeResult1 = applyTimeWindowing( - AscanBlock, sampleRate, distBlock, sosWaterBlock, - expectedSOSWater, offsetElectronicSamples, detectionWindowSOS, + AscanBlock, sampleRate, distBlock, + aSOSWater, offsetElectronicSamples, detectionWindowSOS, minSpeedOfSound, maxSpeedOfSound, gaussWindow); timeResult2 = applyTimeWindowing( - AscanRefBlock, sampleRate, distRefBlock, sosWaterRefBlock, - expectedSOSWater, offsetElectronicSamples, detectionWindowSOS, + AscanRefBlock, sampleRate, distRefBlock, + aSOSWaterRef, offsetElectronicSamples, detectionWindowSOS, minSpeedOfSound, maxSpeedOfSound, gaussWindow); diffStartSearch = timeResult1.startSearch - timeResult2.startSearch; @@ -307,13 +304,13 @@ namespace Recon { if (useTimeWindowing) { shiftInSamples = shiftInSamples - diffStartSearch; } - auto tof = shiftInSamples / sampleRate + distBlock / sosWaterBlock; + auto tof = shiftInSamples / sampleRate + distBlock / aSOSWater; auto sosValue = distBlock / tof; DetectResult result; result.sosValue = sosValue; - auto tofRel = tof - distBlock / sosWaterBlock; + auto tofRel = tof - distBlock / aSOSWater; result.att = detectAttVectorized( - _AscanBlock, _AscanRefBlock, distRefBlock, sosWaterRefBlock, + _AscanBlock, _AscanRefBlock, distRefBlock, aSOSWaterRef, tof, aScanReconstructionFrequency, offsetElectronic, detectionWindowATT); result.tof = tofRel; @@ -423,11 +420,8 @@ namespace Recon { const Aurora::CudaMatrix &AscanRefBlock, const Aurora::CudaMatrix &distBlock, const Aurora::CudaMatrix &distRefBlock, - const Aurora::Matrix &sosWaterBlock, - const Aurora::Matrix &sosWaterRefBlock, - float expectedSOSWater) { - auto _sosWaterBlock = temperatureToSoundSpeed(sosWaterBlock, "marczak").toDeviceMatrix(); - auto _sosWaterRefBlock = temperatureToSoundSpeed(sosWaterRefBlock, "marczak").toDeviceMatrix(); + float aSOSWater, float aSOSWaterRef) + { switch (Recon::transParams::version) { // case 1: { // return detectTofAndAttMex( @@ -442,8 +436,8 @@ namespace Recon { default: auto r = detectTofAndAtt( AscanBlock, AscanRefBlock, distBlock, distRefBlock, - _sosWaterBlock, _sosWaterRefBlock, Recon::transParams::resampleFactor, Recon::transParams::nThreads, - expectedSOSWater, Recon::transParams::useTimeWindowing, + 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); diff --git a/src/transmissionReconstruction/detection/detection.cu b/src/transmissionReconstruction/detection/detection.cu index 8b52570..571c4db 100644 --- a/src/transmissionReconstruction/detection/detection.cu +++ b/src/transmissionReconstruction/detection/detection.cu @@ -66,7 +66,7 @@ CudaMatrix Recon::calculateAttenuationCuda(const CudaMatrix &ascans, CudaMatrix Recon::detectAttVectorizedCuda(const CudaMatrix &Ascan, const CudaMatrix &AscanRef, const CudaMatrix &distRef, - const CudaMatrix &sosWaterRef, + float sosWaterRef, const CudaMatrix &tof, int aScanReconstructionFrequency, float offsetElectronic, int detectionWindowATT) { auto sizeAscan = size(Ascan); @@ -169,19 +169,18 @@ __global__ void guassWindowKernel(float* aStartSearch,float* aEndSearch, } Recon::TimeWindowResultC Recon::applyTimeWindowing(const Aurora::CudaMatrix &AscanBlock, float sampleRate, - const Aurora::CudaMatrix &distBlock, const Aurora::CudaMatrix &sosBlock, - float expectedSOSWater, float startOffset, float segmentLenOffset, + const Aurora::CudaMatrix &distBlock, + float aSOSWater, float startOffset, float segmentLenOffset, float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow) { - auto sosOffset = calculateSOSOffset(sosBlock, expectedSOSWater, distBlock, sampleRate); - + Aurora::CudaMatrix sosOffset = Aurora::zerosCuda(1,1); 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; + auto expectedPosWater = (distBlock / aSOSWater) * sampleRate + startOffset; guassWindowKernel<<>>(calcResult.startSearch.getData(), calcResult.endSearch.getData(), AscanBlock.getData(), AscanBlockProcessed.getData(), expectedPosWater.getData(), AscanBlock.getDimSize(0)); @@ -253,9 +252,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 expectedSOSWater, + int resampleFactor,int nthreads, float aSOSWater, float aSOSWaterRef, int useTimeWindowing, int aScanReconstructionFrequency,int detectionWindowATT, float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow) @@ -269,12 +266,12 @@ Recon::DetectResultC Recon::detectTofAndAtt( timeResult2.AscanBlockProcessed = AscanRefBlock; if (useTimeWindowing == 1) { timeResult1 = applyTimeWindowing( - AscanBlock, sampleRate, distBlock, sosWaterBlock, - expectedSOSWater, offsetElectronicSamples, detectionWindowSOS, + AscanBlock, sampleRate, distBlock, + aSOSWater, offsetElectronicSamples, detectionWindowSOS, minSpeedOfSound, maxSpeedOfSound, gaussWindow); timeResult2 = applyTimeWindowing( - AscanRefBlock, sampleRate, distRefBlock, sosWaterRefBlock, - expectedSOSWater, offsetElectronicSamples, detectionWindowSOS, + AscanRefBlock, sampleRate, distRefBlock, + aSOSWaterRef, offsetElectronicSamples, detectionWindowSOS, minSpeedOfSound, maxSpeedOfSound, gaussWindow); diffStartSearch = timeResult1.startSearch - timeResult2.startSearch; @@ -296,8 +293,8 @@ Recon::DetectResultC Recon::detectTofAndAtt( CudaMatrix c_1_1; { auto x = fft(_AscanBlock, m2); - auto y = fft(_AscanRefBlock, m2); - c_1_1 = x * conj(y); + auto y = conj(fft(_AscanRefBlock, m2)); + c_1_1 = x * y; } c_1_2 = ifft(c_1_1); } @@ -313,13 +310,13 @@ Recon::DetectResultC Recon::detectTofAndAtt( if (useTimeWindowing) { shiftInSamples = shiftInSamples - diffStartSearch; } - auto tof = shiftInSamples / sampleRate + distBlock / sosWaterBlock; + auto tof = shiftInSamples / sampleRate + distBlock / aSOSWater; auto sosValue = distBlock / tof; Recon::DetectResultC result; result.sosValue = sosValue; - auto tofRel = tof - distBlock / sosWaterBlock; + auto tofRel = tof - distBlock / aSOSWater; result.att = detectAttVectorizedCuda( - _AscanBlock, _AscanRefBlock, distRefBlock, sosWaterRefBlock, + _AscanBlock, _AscanRefBlock, distRefBlock, aSOSWaterRef, tof, aScanReconstructionFrequency, offsetElectronic, detectionWindowATT); result.tof = tofRel; diff --git a/src/transmissionReconstruction/detection/detection.cuh b/src/transmissionReconstruction/detection/detection.cuh index 0767e9e..75e841d 100644 --- a/src/transmissionReconstruction/detection/detection.cuh +++ b/src/transmissionReconstruction/detection/detection.cuh @@ -25,13 +25,13 @@ struct TimeWindowResultC { CudaMatrix detectAttVectorizedCuda(const CudaMatrix &Ascan, const CudaMatrix &AscanRef, const CudaMatrix &distRef, - const CudaMatrix &sosWaterRef, + float 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, + const Aurora::CudaMatrix &distBlock, + float aSOSWater, float startOffset, float segmentLenOffset, float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow); SearchPositionC calculateStarEndSearchPosition(const CudaMatrix &aVDistBlock, float minSpeedOfSound, float maxSpeedOfSound, @@ -42,9 +42,7 @@ struct TimeWindowResultC { 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 resampleFactor,int nthreads, float aSOSWater, float aSOSWaterRef, int useTimeWindowing, int aScanReconstructionFrequency,int detectionWindowATT, float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow); diff --git a/src/transmissionReconstruction/detection/detection.h b/src/transmissionReconstruction/detection/detection.h index 6f104b7..359d05a 100644 --- a/src/transmissionReconstruction/detection/detection.h +++ b/src/transmissionReconstruction/detection/detection.h @@ -32,14 +32,14 @@ calculateStarEndSearchPosition(const Aurora::Matrix &aVDistBlock, TimeWindowResult applyTimeWindowing( const Aurora::Matrix &AscanBlock, float sampleRate, - const Aurora::Matrix &distBlock, const Aurora::Matrix &sosBlock, - float expectedSOSWater, float startOffset, float segmentLenOffset, + const Aurora::Matrix &distBlock, float sosWater, + float startOffset, float segmentLenOffset, float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow); Aurora::Matrix detectAttVectorized(const Aurora::Matrix &Ascan, const Aurora::Matrix &AscanRef, const Aurora::Matrix &distRef, - const Aurora::Matrix &sosWaterRef, + float sosWaterRef, const Aurora::Matrix &tof, int aScanReconstructionFrequency, float offsetElectronic, int detectionWindowATT); @@ -47,8 +47,7 @@ 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 expectedSOSWater, + float aSOSWater, float aSOSWaterRef, int useTimeWindowing, int aScanReconstructionFrequency, float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow); @@ -57,8 +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 expectedSOSWater, + int resampleFactor, int nthreads, float aSOSWater, float aSOSWaterRef, int useTimeWindowing, int aScanReconstructionFrequency, int detectionWindowATT, float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow); @@ -79,7 +77,7 @@ DetectResult transmissionDetection( const Aurora::CudaMatrix &AscanBlock, const Aurora::CudaMatrix &AscanRefBlock, const Aurora::CudaMatrix &distBlock, const Aurora::CudaMatrix &distRefBlock, - const Aurora::Matrix &sosWaterBlock, const Aurora::Matrix &sosWaterRefBlock, float expectedSOSWater); + float aSOSWater, float aSOSWaterRef); } // namespace Recon diff --git a/src/transmissionReconstruction/detection/getTransmissionData.cpp b/src/transmissionReconstruction/detection/getTransmissionData.cpp index 3b993e1..62782f6 100644 --- a/src/transmissionReconstruction/detection/getTransmissionData.cpp +++ b/src/transmissionReconstruction/detection/getTransmissionData.cpp @@ -1,4 +1,5 @@ #include "getTransmissionData.h" +#include "common/getMeasurementMetaData.h" #include "getTransmissionData.cuh" #include "AuroraDefs.h" #include "CudaMatrix.h" @@ -45,14 +46,12 @@ namespace Matrix attData; Matrix senderBlock; Matrix receiverBlock; - Matrix waterTempBlock; MetaInfos metaInfos; CudaMatrix ascanBlock; CudaMatrix ascanBlockRef; Matrix dists; Matrix distRefBlock; - Matrix waterTempRefBlock; }; std::map BLOCK_OF_TRANSIMISSIONDARA_BUFFER; @@ -79,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 TasTemps& aTasTemps, const Matrix& aExpectedSOSWater, GeometryInfo aGeom, GeometryInfo& aGeomRef, + const TempInfo& aTemp, const TempInfo& aTempRef, GeometryInfo aGeom, GeometryInfo& aGeomRef, const Matrix& aSnrRmsNoise, const Matrix& aSnrRmsNoiseRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef, const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef) { @@ -135,9 +134,6 @@ namespace Matrix dists = distanceBetweenTwoPoints(blockData.senderPositionBlock, blockData.receiverPositionBlock); Matrix distRefBlock = distanceBetweenTwoPoints(blockDataRef.senderPositionBlock, blockDataRef.receiverPositionBlock); - - Matrix waterTempBlock = calculateWaterTemperature(aTasTemps.waterTempPreCalc_sl, aTasTemps.waterTempPreCalc_rl, blockData.slBlock, blockData.rlBlock, blockData.mpBlock); - Matrix waterTempRefBlock = calculateWaterTemperature(aTasTemps.waterTempRefPreCalc_sl, aTasTemps.waterTempRefPreCalc_rl, blockData.slBlock, blockData.rlBlock, blockDataRef.mpBlock); if(transParams::saveDetection || transParams::outlierOnTasDetection || transParams::saveDebugInfomation) { @@ -150,13 +146,11 @@ namespace result.metaInfos = metaInfos; result.senderBlock = blockData.senderPositionBlock; result.receiverBlock = blockData.receiverPositionBlock; - result.waterTempBlock = waterTempBlock; result.ascanBlock = ascanBlock; result.ascanBlockRef = ascanBlockRef; result.dists = dists; result.distRefBlock = distRefBlock; - result.waterTempRefBlock = waterTempRefBlock; // DetectResult detect = transmissionDetection(ascanBlock, ascanBlockRef, dists, distRefBlock, waterTempBlock, waterTempRefBlock, aExpectedSOSWater[0]); // result.attData = detect.att; @@ -168,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 TasTemps& aTasTemps, const Matrix& aExpectedSOSWater, GeometryInfo aGeom, GeometryInfo aGeomRef, + const TempInfo& aTemp, const TempInfo& aTempRef, 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, aTasTemps, - aExpectedSOSWater, aGeom, aGeomRef, aSnrRmsNoise, aSnrRmsNoiseRef, + auto buffer = getBlockOfTransmissionData(aMp, aMpRef, aSl, aSn, aRlList, aRnList, aTemp, + aTempRef, aGeom, aGeomRef, aSnrRmsNoise, aSnrRmsNoiseRef, aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef); std::unique_lock lock(CREATE_BUFFER_MUTEX); BLOCK_OF_TRANSIMISSIONDARA_BUFFER[std::to_string(aIndex)] = buffer; @@ -184,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 TasTemps& aTasTemps, const Matrix& aExpectedSOSWater, GeometryInfo aGeom, GeometryInfo aGeomRef, + const TempInfo& aTemp, const TempInfo& aTempRef, GeometryInfo aGeom, GeometryInfo aGeomRef, const Matrix& aSnrRmsNoise, const Matrix& aSnrRmsNoiseRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef, const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef) { @@ -208,7 +202,7 @@ void createThreadForGetBlockOfTransmissionData(const Matrix& aMotorPos, const M CREATE_BUFFER_CONDITION.wait(lock, []{return BUFFER_COUNT lockBufferCount(CREATE_BUFFER_MUTEX); @@ -353,7 +339,6 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con receiverList = removeDataFromArrays(receiverList, filter); tofDataTotal = removeDataFromArrays(tofDataTotal, filter); attDataTotal = removeDataFromArrays(attDataTotal, filter); - waterTempList = removeDataFromArrays(waterTempList, filter); Matrix valid; if(transParams::applyCalib) @@ -383,7 +368,6 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con attDataTotal = removeDataFromArrays(attDataTotal, valid); senderList = removeDataFromArrays(senderList, valid); receiverList = removeDataFromArrays(receiverList, valid); - waterTempList = removeDataFromArrays(waterTempList, valid); dataInfno.numPossibleScans = numData; dataInfno.numValidScans = sum(valid); @@ -439,6 +423,5 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con result.receiverList = receiverList; result.senderList = senderList; result.dataInfo = dataInfno; - result.waterTempList = waterTempList; return result; } \ No newline at end of file diff --git a/src/transmissionReconstruction/detection/getTransmissionData.h b/src/transmissionReconstruction/detection/getTransmissionData.h index 46c3a1d..431e733 100644 --- a/src/transmissionReconstruction/detection/getTransmissionData.h +++ b/src/transmissionReconstruction/detection/getTransmissionData.h @@ -37,7 +37,6 @@ namespace Recon Aurora::Matrix attDataTotal; Aurora::Matrix senderList; Aurora::Matrix receiverList; - Aurora::Matrix waterTempList; DataInfo dataInfo; }; diff --git a/src/transmissionReconstruction/startTransmissionReconstruction.cpp b/src/transmissionReconstruction/startTransmissionReconstruction.cpp index 688b4bf..ec804db 100644 --- a/src/transmissionReconstruction/startTransmissionReconstruction.cpp +++ b/src/transmissionReconstruction/startTransmissionReconstruction.cpp @@ -3,6 +3,7 @@ #include "Function2D.h" #include "MatlabWriter.h" +#include "Matrix.h" #include "detection/getTransmissionData.h" #include "log/log.h" #include "config/config.h" @@ -26,9 +27,8 @@ TransmissionReconstructionResult Recon::startTransmissionReconstruction(const Au auto transmissionData = getTransmissionData(aMotorPos, aMotoPosRef, aSlList, aSnList, aRlList, aRnList, aTemp, aTempRef, aGeom, aGeomRef, aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef); Matrix dists = Recon::distanceBetweenTwoPoints(transmissionData.senderList, transmissionData.receiverList); - Matrix sosRef = Recon::temperatureToSoundSpeed(transmissionData.waterTempList, "marczak"); - Matrix valid = Recon::checkTofDetections(transmissionData.tofDataTotal, dists, sosRef, + Matrix valid = Recon::checkTofDetections(transmissionData.tofDataTotal, dists, aTemp.expectedSOSWater[0], Recon::transParams::minSpeedOfSound,Recon::transParams::maxSpeedOfSound).valid; if(transParams::qualityCheck) { diff --git a/test/DataFilter_Test.cpp b/test/DataFilter_Test.cpp index b9b0468..9920173 100644 --- a/test/DataFilter_Test.cpp +++ b/test/DataFilter_Test.cpp @@ -74,28 +74,28 @@ TEST_F(DataFilter_Test, filterTransmissionAngle) { } -TEST_F(DataFilter_Test, checkTofDetections) { - MatlabReader m("/home/sun/testData/checkTofDetections.mat"); - auto receiverList = m.read("receiverList"); - auto senderList = m.read("senderList"); - auto tofDataTotal = m.read("tofDataTotal"); - auto waterTempList = m.read("waterTempList"); - auto tofValues = m.read("tofValues"); - auto valid = m.read("valid"); - Aurora::Matrix dists = Recon::distanceBetweenTwoPoints(senderList, receiverList); - Aurora::Matrix sosRef = Recon::temperatureToSoundSpeed(waterTempList, "marczak"); - auto result = Recon::checkTofDetections(tofDataTotal, dists, sosRef, Recon::transParams::minSpeedOfSound,Recon::transParams::maxSpeedOfSound); +// TEST_F(DataFilter_Test, checkTofDetections) { +// MatlabReader m("/home/sun/testData/checkTofDetections.mat"); +// auto receiverList = m.read("receiverList"); +// auto senderList = m.read("senderList"); +// auto tofDataTotal = m.read("tofDataTotal"); +// auto waterTempList = m.read("waterTempList"); +// auto tofValues = m.read("tofValues"); +// auto valid = m.read("valid"); +// Aurora::Matrix dists = Recon::distanceBetweenTwoPoints(senderList, receiverList); +// Aurora::Matrix sosRef = Recon::temperatureToSoundSpeed(waterTempList, "marczak"); +// auto result = Recon::checkTofDetections(tofDataTotal, dists, sosRef, Recon::transParams::minSpeedOfSound,Recon::transParams::maxSpeedOfSound); - for (size_t i = 0; i < result.valid.getDataSize(); i++) - { - EXPECT_FLOAT_AE(valid.getData()[i],result.valid.getData()[i]) << " :"<