diff --git a/src/common/dataBlockCreation/blockingGeometryInfo.cpp b/src/common/dataBlockCreation/blockingGeometryInfo.cpp index afee939..c3e58f0 100644 --- a/src/common/dataBlockCreation/blockingGeometryInfo.cpp +++ b/src/common/dataBlockCreation/blockingGeometryInfo.cpp @@ -25,6 +25,7 @@ GeometryBlock Recon::blockingGeometryInfos(const GeometryInfo& aGeom, const Auro sizeData[3] = aGeom.receiverPositions.size(); size_t matrixSize = aGeom.receiverPositions[0].getDataSize(); Matrix receiverPositionsSize = Matrix::New(sizeData, 1, 4); + #pragma omp parallel for for(int i=0; i 0) + { + result.ascanBlockPreprocessed = preprocessAscanBlock(ascanBlock.ascanBlock, aMeasInfo); + } + else + { + result.ascanBlockPreprocessed = ascanBlock.ascanBlock; + } + + return result; +} \ No newline at end of file diff --git a/src/common/dataBlockCreation/getAScanBlockPreprocessed.h b/src/common/dataBlockCreation/getAScanBlockPreprocessed.h new file mode 100644 index 0000000..7285e1c --- /dev/null +++ b/src/common/dataBlockCreation/getAScanBlockPreprocessed.h @@ -0,0 +1,31 @@ +#ifndef GETASCANBLOCK_PREPROCESSED_H +#define GETASCANBLOCK_PREPROCESSED_H + +#include "Matrix.h" +#include "src/common/getGeometryInfo.h" +#include "src/common/getMeasurementMetaData.h" + +class Parser; + +namespace Recon +{ + struct AscanBlockPreprocessed + { + Aurora::Matrix ascanBlockPreprocessed; + Aurora::Matrix mpBlock; + Aurora::Matrix slBlock; + Aurora::Matrix snBlock; + Aurora::Matrix rlBlock; + Aurora::Matrix rnBlock; + Aurora::Matrix senderPositionBlock; + Aurora::Matrix receiverPositionBlock; + Aurora::Matrix gainBlock; + }; + + AscanBlockPreprocessed getAscanBlockPreprocessed(Parser* aParser, const Aurora::Matrix& aMp, const Aurora::Matrix& aSl, const Aurora::Matrix& aSn, + const Aurora::Matrix& aRl, const Aurora::Matrix& aRn, GeometryInfo& aGeom, const MeasurementInfo& aMeasInfo, + bool aApplyFilter, bool aTransReco); +} + + +#endif \ No newline at end of file diff --git a/src/common/dataBlockCreation/getAscanBlock.cpp b/src/common/dataBlockCreation/getAscanBlock.cpp index 2f239f9..f8c5ab7 100644 --- a/src/common/dataBlockCreation/getAscanBlock.cpp +++ b/src/common/dataBlockCreation/getAscanBlock.cpp @@ -11,6 +11,8 @@ #include "Data/ElectricIndex.h" #include "Data/MetaData.h" #include "ShotList/ShotList.h" +#include +#include #include using namespace Recon; @@ -18,12 +20,17 @@ using namespace Aurora; namespace { - std::vector findTasIndex(TasIndicesPointer aTasIndices, int aTasNum) + std::vector findTasAndElementIndex(TasIndicesPointer aTasIndices, ReceiverIndicesPointer aReceiverIndices, const Matrix& aRl, const Matrix& aRn) { std::vector result; + Matrix sortRl = Matrix::copyFromRawData(aRl.getData(), aRl.getDimSize(0), aRl.getDimSize(1), aRl.getDimSize(2)); + Matrix sortRn = Matrix::copyFromRawData(aRn.getData(), aRn.getDimSize(0), aRn.getDimSize(1), aRn.getDimSize(2)); + std::sort(sortRl.getData(), sortRl.getData() + sortRl.getDataSize()); + std::sort(sortRn.getData(), sortRn.getData() + sortRn.getDataSize()); for(int i=0; igetMetaData().getSampleNumber(); double* ascanBlockData = new double[ascanBlockSize]; double* mpBlockData = new double[numScans]; @@ -57,30 +64,39 @@ AscanBlock Recon::getAscanBlock(Parser* aParser, const Aurora::Matrix& aMp, cons OneTasAScanData oenTasData = aParser->getOneTasAscanDataOfMotorPosition(1); auto tasIndices = aParser->getMetaData().getTasIndices(); auto receiverIndices = aParser->getMetaData().getReceiverIndices(); + auto tasElementMapper = findTasAndElementIndex(tasIndices, receiverIndices, aRl, aRn); + for(int mpIndex=0; mpIndexgetOneTasAscanDataOfMotorPosition(aSl[slIndex]); + for(int snIndex=0; snIndexsearchAscanDataFromOneTasAscanDataOfMP(oenTasData, ElementIndex(GeometryIndex(aSn[snIndex])), TasElementIndex(aRl[rlIndex], ElementIndex(GeometryIndex(receiverIndices.get()[rElementIndex[rnIndex]]))), mp); - std::copy(ascan.get() ,ascan.get() + ascan.getAscanDataLength(), ascanBlockData); - ascanBlockData += ascan.getAscanDataLength(); - + AScanData ascan = aParser->searchAscanDataFromOneTasAscanDataOfMP(oenTasData, ElementIndex(GeometryIndex(aSn[snIndex])), TasElementIndex(tasIndices.get()[tasElementMapper[mapperIndex]], ElementIndex(GeometryIndex(receiverIndices.get()[tasElementMapper[mapperIndex]]))), mp); + double* startPointer = ascanBlockData + numScansIndex * ascan.getAscanDataLength(); + std::copy(ascan.get() ,ascan.get() + ascan.getAscanDataLength(), startPointer); + //ascanBlockData += ascan.getAscanDataLength(); mpBlockData[numScansIndex] = aMp[mpIndex]; slBlockData[numScansIndex] = aSl[slIndex]; snBlockData[numScansIndex] = aSn[snIndex]; - rlBlockData[numScansIndex] = aRl[rlIndex]; - rnBlockData[numScansIndex] = receiverIndices.get()[rElementIndex[rnIndex]]; + rlBlockData[numScansIndex] = tasIndices.get()[tasElementMapper[mapperIndex]]; + rnBlockData[numScansIndex] = receiverIndices.get()[tasElementMapper[mapperIndex]]; gainBlockData[numScansIndex] = ascan.getAmplification()[0]; - ++numScansIndex; + //++numScansIndex; + //++mapperIndex; } } } diff --git a/src/common/precalculateChannelList.cpp b/src/common/precalculateChannelList.cpp new file mode 100644 index 0000000..f8c04dd --- /dev/null +++ b/src/common/precalculateChannelList.cpp @@ -0,0 +1,49 @@ +#include "precalculateChannelList.h" +#include "Function.h" +#include "Matrix.h" + +using namespace Aurora; +using namespace Recon; + +namespace +{ + double USCTTAS2DAQChannels(double aRl, double aRn, const MeasurementInfo& aExpInfo, const PreComputes& aPreComputes) + { + size_t size = aPreComputes.measuredCE_receiverIndices.getDataSize(); + if(aExpInfo.Hardware == "USCT3dv2") + { + //暂不实现 + } + else if(aExpInfo.Hardware == "USCT3dv3") + { + if (aPreComputes.measuredCEUsed == 1) + { + for(int i=0; i + +Aurora::Matrix Recon::preprocessAscanBlock(const Aurora::Matrix& aAscans, const MeasurementInfo& aMeasInfo) +{ + Aurora::Matrix result = aAscans; + size_t size = aAscans.getDataSize(); + short* ascanData = new short[size]; + std::copy(aAscans.getData(), aAscans.getData() + size, ascanData); + if(aMeasInfo.ascanDataType == "float16") + { + result = Aurora::convertfp16tofloat(ascanData, aAscans.getDimSize(0), aAscans.getDimSize(1)); + } + + //暂不考虑实现二代逻辑 + // if isfield(measInfo, 'Bandpassundersampling') && (measInfo.Bandpassundersampling == 1) + // AScans = reconstructBandpasssubsampling(AScans, params.aScanReconstructionFrequency, measInfo.SampleRate); + // end + + return result; +} \ No newline at end of file diff --git a/src/common/preprocessAscanBlock.h b/src/common/preprocessAscanBlock.h new file mode 100644 index 0000000..961ffa0 --- /dev/null +++ b/src/common/preprocessAscanBlock.h @@ -0,0 +1,12 @@ +#ifndef PREPROCESS_ASCANBLOCK_H +#define PREPROCESS_ASCANBLOCK_H + +#include "Matrix.h" +#include "getMeasurementMetaData.h" + +namespace Recon +{ + Aurora::Matrix preprocessAscanBlock(const Aurora::Matrix& aAscans, const MeasurementInfo& aMeasInfo); +} + +#endif \ No newline at end of file diff --git a/src/common/qualityReview.cpp b/src/common/qualityReview.cpp new file mode 100644 index 0000000..334cfe8 --- /dev/null +++ b/src/common/qualityReview.cpp @@ -0,0 +1,11 @@ +#include "qualityReview.h" +#include "Matrix.h" + +using namespace Recon; +using namespace Aurora; + +void Recon::qualityReview(double aNumValiData, double aNumTotalData) +{ + double score = aNumValiData / aNumTotalData; + //to do, 输出报警用 +} \ No newline at end of file diff --git a/src/common/qualityReview.h b/src/common/qualityReview.h new file mode 100644 index 0000000..6a14a6a --- /dev/null +++ b/src/common/qualityReview.h @@ -0,0 +1,11 @@ +#ifndef QUALITY_REVIEW_H +#define QUALITY_REVIEW_H + +#include "Matrix.h" + +namespace Recon +{ + void qualityReview(double aNumValiData, double aNumTotalData); +} + +#endif \ No newline at end of file diff --git a/src/common/temperatureCalculation/extractTasTemperature.cpp b/src/common/temperatureCalculation/extractTasTemperature.cpp index 95f4736..21a3df9 100644 --- a/src/common/temperatureCalculation/extractTasTemperature.cpp +++ b/src/common/temperatureCalculation/extractTasTemperature.cpp @@ -25,7 +25,7 @@ namespace { if(aReferenceTemps.getDimSize(1) >= aMotorPos) { - result = mean(aReferenceTemps($,aMotorPos).toMatrix(), FunctionDirection::All, false)[0]; + result = mean(aReferenceTemps($,aMotorPos - 1).toMatrix(), FunctionDirection::All, false)[0]; } else { @@ -61,6 +61,5 @@ Matrix Recon::extractTasTemperature(const Matrix& aTasTemperature, const Matrix& } } Matrix result = Matrix::New(resultData, aTasList.getDataSize(), aMotorPosList.getDataSize()); - result.printf(); return result; } \ No newline at end of file diff --git a/src/startReconstructions.cpp b/src/startReconstructions.cpp new file mode 100644 index 0000000..a984dbd --- /dev/null +++ b/src/startReconstructions.cpp @@ -0,0 +1,215 @@ +#include "startReconstructions.h" +#include "Data/TemperatureData.h" +#include "Function.h" +#include "config/config.h" +#include "common/getMeasurementMetaData.h" +#include "common/getGeometryInfo.h" +#include "common/estimatePulseLength.h" +#include "transmissionReconstruction/startTransmissionReconstruction.h" + +#include +#include +#include "Parser/Parser.h" +#include "Parser/Data/MetaData.h" + +#include "Matrix.h" +#include "Function1D.h" +#include "Function2D.h" +#include "src/common/ceMatchedFilterHandling.h" +#include "src/reflectionReconstruction/preprocessData/estimateOffset.h" + + +using namespace Recon; +using namespace Aurora; + +void Recon::startReconstructions() +{ + // std::string dataPath = "/home/AScans_Data/volunteer_20230328/20230328T123237/"; + // std::string refPath = "/home/AScans_Data/volunteer_20230328/20230328T123237/"; + std::string dataPath = "/home/AScans_Data/ADW_TAS_Issue/20230512T135108/"; + std::string refPath = "/home/AScans_Data/ADW_TAS_Issue/20230512T141626/"; + std::string outputPath; + Parser dataParser(dataPath); + Parser refParser(refPath); + std::string rootMeasUniqueID = dataParser.getMetaData().getMeasurementID(); + std::string rootRefUniqueID = refParser.getMetaData().getMeasurementID(); + + //init total used receiver/emitter/motorPos + Matrix motorPosTotal, slList, snList, rlList, rnList; + if(transParams::runTransmissionReco && reflectParams::runReflectionReco) + { + motorPosTotal = auroraUnion(reflectParams::motorPos, transParams::motorPos); + slList = auroraUnion(reflectParams::senderTasList, transParams::senderTasList); + snList = auroraUnion(reflectParams::senderElementList, transParams::senderElementList); + rlList = auroraUnion(reflectParams::receiverTasList, transParams::receiverTasList); + rnList = auroraUnion(reflectParams::receiverElementList, transParams::receiverElementList); + } + else if(reflectParams::runReflectionReco) + { + motorPosTotal = reflectParams::motorPos; + slList = reflectParams::senderTasList; + snList = reflectParams::senderElementList; + rlList = reflectParams::receiverTasList; + rnList = reflectParams::receiverElementList; + } + else if(transParams::runTransmissionReco) + { + motorPosTotal = transParams::motorPos; + slList = transParams::senderTasList; + snList = transParams::senderElementList; + rlList = transParams::receiverTasList; + rnList = transParams::receiverElementList; + } + else + { + return; + } + + //getMeasurementMetaData + double maxNumTAS = max(auroraUnion(slList, rlList)).getData()[0]; + MeasurementInfo expInfo = loadMeasurementInfos(&dataParser); + TempInfo temp = getTemperatureInfo(&dataParser, maxNumTAS); + CEInfo ce = getCEInfo(&dataParser, expInfo); + 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); + tempRef = getTemperatureInfo(&refParser, maxNumTAS); + ceRef = getCEInfo(&refParser, expInfoRef); + transformationInfo = getTransformationMatrix(&refParser, motorPosTotal); + transformationMatricesRef = transformationInfo.rotationMatrix; + motorPosAvailableRef = transformationInfo.motorPos; + if(transformationMatricesRef.isNull()) + { + Matrix motorPos1 = Matrix::fromRawData(new double[1] {1}, 1); + transformationMatricesRef = getTransformationMatrix(&refParser, motorPos1).rotationMatrix; + } + else + { + Matrix mpRef_inter = intersect(motorPosAvailableRef, transParams::motorPos); + //extractMatchingReferenceMotorPosition + for(int i=0; i +#include + +#include "MatlabReader.h" +#include +#include +#include using namespace Recon; using namespace Aurora; -void Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, const Aurora::Matrix& aMotoPosRef, const Aurora::Matrix& aSlList, +namespace +{ + struct BlockOfTransmissionData + { + Matrix tofData; + Matrix attData; + Matrix senderBlock; + Matrix receiverBlock; + Matrix waterTempBlock; + MetaInfos metaInfos; + }; + + Matrix prepareAScansForTransmissionDetection(const Matrix& aAscanBlock, const Matrix& aGainBlock) + { + Matrix result = aAscanBlock / repmat(aGainBlock, aAscanBlock.getDimSize(0), 1); + result = result - repmat(mean(result,FunctionDirection::Column), result.getDimSize(0), 1); + return result; + } + + 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 Matrix& aSnrRmsNoise, const Matrix& aSnrRmsNoiseRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef, + const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef) + { + BlockOfTransmissionData result; + MetaInfos metaInfos; + auto blockData = getAscanBlockPreprocessed(aParser, aMp, aSl, aSn, aRlList, aRnList, aGeom, aExpInfo, true, true); + auto blockDataRef = getAscanBlockPreprocessed(aParserRef, aMpRef, aSl, aSn, aRlList, aRnList, aGeomRef, aExpInfoRef, true, true); + Matrix ascanBlock = prepareAScansForTransmissionDetection(blockData.ascanBlockPreprocessed, blockData.gainBlock); + Matrix ascanBlockRef = prepareAScansForTransmissionDetection(blockDataRef.ascanBlockPreprocessed, blockDataRef.gainBlock); + if(aExpInfo.Hardware == "USCT3dv3") + { + Matrix channelList = precalculateChannelList(aRlList, aRnList, aExpInfo, aPreComputes); + double* channelListSizeData = Aurora::malloc(2); + channelListSizeData[0] = channelList.getDimSize(0); + channelListSizeData[1] = channelList.getDimSize(1); + Matrix channelListSize = Matrix::New(channelListSizeData, 2, 1); + Matrix ind = sub2ind(channelListSize, {blockData.rlBlock, blockData.rnBlock}); + size_t channelListBlockSize = ind.getDataSize(); + double* channelListBlockData = Aurora::malloc(channelListBlockSize); + for(size_t i=0; i 2 % i.e. detection folder exists and is not empty // % Load transmission detection data @@ -30,4 +198,204 @@ void Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, const Aurora::M aGeom.sensData = precalcSensitivity(aSlList, aSnList, aRlList, aRnList, aMotorPos, aGeom); aGeomRef.sensData = aGeom.sensData; + Matrix rmsNoise, rmsNoiseRef; + if (transParams::applyCalib) + { + rmsNoise = estimateNoiseValueFromAScans(aSnList, aRnList, aGeom, aExpInfo, aParser); + rmsNoiseRef = estimateNoiseValueFromAScans(aSnList, aRnList, aGeom, aExpInfo, aParserRef); + } + + size_t numScans = aMotorPos.getDataSize() * aSlList.getDataSize() * aSnList.getDataSize() * aRlList.getDataSize() * aRnList.getDataSize(); + Matrix tofDataTotal = Matrix::fromRawData(new double[numScans], 1, numScans) + NAN; + Matrix attDataTotal = Matrix::fromRawData(new double[numScans], 1, numScans) + NAN; + Matrix waterTempList = zeros(1,numScans,1); + Matrix senderList = zeros(3,numScans,1); + Matrix receiverList = zeros(3,numScans,1); + Matrix snrValues = zeros(1,numScans,1); + Matrix snrValuesRef = zeros(1,numScans,1); + + Matrix mpBlockTotal; + Matrix slBlockTotal; + Matrix snBlockTotal; + Matrix rlBlockTotal; + Matrix rnBlockTotal; + if(transParams::saveDetection || transParams::outlierOnTasDetection || transParams::saveDebugInfomation) + { + mpBlockTotal = zeros(1,numScans,1); + slBlockTotal = zeros(1,numScans,1); + snBlockTotal = zeros(1,numScans,1); + rlBlockTotal = zeros(1,numScans,1); + rnBlockTotal = zeros(1,numScans,1); + } + + int numData = 0; + int numPossibleScans = 0; + + for(int i=0; i