add getMeasurementMetaData

src/common/getMeasurementMetaData.cpp

@@ -0,0 +1,335 @@
+#include "getMeasurementMetaData.h"
+#include "Matrix.h"
+#include "ceMatchedFilterHandling.h"
+
+#include "ShotList/ShotList.h"
+#include "Parser.h"
+#include "Data/MetaData.h"
+#include "Data/MovementData.h"
+#include "Data/TemperatureData.h"
+#include "Data/CEMeasuredData.h"
+#include "Data/CEData.h"
+
+#include "Function.h"
+#include "Function1D.h"
+#include "Function2D.h"
+#include "AuroraDefs.h"
+#include "mkl.h"
+
+#include "../config/config.h"
+#include <algorithm>
+
+using namespace Recon;
+using namespace Aurora;
+
+Matrix Recon::getAvailableMotorPositions(Parser* aParser)
+{
+    std::vector<double> mpData;
+    if(aParser->getShotList()->hasCEMeasured())
+    {
+        mpData.push_back(0);
+    }
+    int motorPosition = aParser->getMetaData().getAperturePositionNumber();
+    for(int i=1; i<=motorPosition; ++i)
+    {
+        mpData.push_back(i);
+    }
+
+    return Matrix::copyFromRawData(mpData.data(), mpData.size());
+}
+
+//已验证 完全正确
+MeasurementInfo Recon::loadMeasurementInfos(Parser* aParser)
+{
+    MeasurementInfo result;
+    MetaData metaData = aParser->getMetaData();
+    result.numberSamples = metaData.getSampleNumber();
+    result.Bandpassundersampling = false;
+    result.sampleRate = metaData.getSamplerate();
+    result.ascanDataType = metaData.getDataType();
+    result.Hardware = "USCT3dv3";
+    result.HardwareVersion = "3.0";
+    result.EOffset = 0;
+    result.Wavelength = 3;
+    result.expectedAScanLength = result.numberSamples;
+    result.rootMeasUniqueID = metaData.getMeasurementID();
+    if(result.Bandpassundersampling)
+    {
+        result.expectedAScanLength = reflectParams::aScanReconstructionFrequency / result.sampleRate * result.expectedAScanLength;
+    }
+    return result;
+}
+
+//已验证 完全正确
+TransFormInfo Recon::getTransformationMatrix(Parser* aParser, const Matrix& aMotorPosList)
+{
+    TransFormInfo result;
+    bool movementRealAvailable = true;
+    Matrix transformationMatrixList;
+    MovementsListRealPointer listRealData = aParser->getMovementData().getMovementsListReal();
+    int columns = 2;//一个角度一个高度
+    unsigned long long length = listRealData.getLength();
+    int rows = length / columns;
+    double* data = new double[length];
+    std::copy(listRealData.get(), listRealData.get() + length, data);
+    std::vector<int> info = {rows, columns, 1};
+    Matrix movementsListReal(std::shared_ptr<double>(data, std::default_delete<double[]>()), info);
+    double* motorPosListAvailableData = Aurora::malloc(rows);
+    for(int i=0; i<rows; ++i)
+    {
+        motorPosListAvailableData[i] = i + 1;
+    }
+    Matrix motorPosListAvailable = Matrix::New(motorPosListAvailableData, rows);
+    Matrix motorPosList = intersect(motorPosListAvailable, aMotorPosList);
+    result.motorPos = motorPosList;
+
+    if(movementsListReal.isNan())
+    {
+        movementRealAvailable = false;
+    }
+    if(movementRealAvailable)
+    {
+        int maxMotorPosNum = motorPosList.getDimSize(0);
+        double* transformationMatrixListData = Aurora::malloc(16*maxMotorPosNum);
+        double zero = 0.0;
+        cblas_dcopy(16*maxMotorPosNum,&zero,0,transformationMatrixListData,1);
+        for(int i=0; i<maxMotorPosNum; ++i)
+        {
+            double rotation = movementsListReal.getData()[i];
+            double distance = movementsListReal.getData()[i + 2];
+            transformationMatrixListData[16*i + 0] = cos(rotation*PI/180);
+            transformationMatrixListData[16*i + 1] = sin(rotation*PI/180);
+            transformationMatrixListData[16*i + 4] = -sin(rotation*PI/180);
+            transformationMatrixListData[16*i + 5] = cos(rotation*PI/180);
+            transformationMatrixListData[16*i + 10] = 1;
+            transformationMatrixListData[16*i + 14] = distance;
+            transformationMatrixListData[16*i + 15] = 1;
+        }
+        transformationMatrixList = Matrix::New(transformationMatrixListData, 4, 4,maxMotorPosNum);
+    }
+    else
+    {
+        auto rotationMatrix = aParser->getMovementData().getRotationMatrix();
+        size_t maxMotorPosNum = rotationMatrix.size();
+        double* transformationMatrixListData = new double[16*maxMotorPosNum];
+        for(size_t i=0; i<maxMotorPosNum; ++i)
+        {
+            std::copy(rotationMatrix.at(i).get(), rotationMatrix.at(i).get()+16, transformationMatrixListData + i*16);
+        }
+        transformationMatrixList = Matrix::fromRawData(transformationMatrixListData, 4, 4, maxMotorPosNum);
+    }
+    result.rotationMatrix = transformationMatrixList;
+
+    return result;
+}
+
+//已验证，完全正确
+TempInfo Recon::getTemperatureInfo(Parser* aParser, double aNumTas)
+{
+    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();
+    }
+
+    double* jumoTempData = Aurora::malloc(jumoTempNum * 4);
+    for(int i=0; i<jumoTempNum; ++i)
+    {
+        jumoTempData[4*i] = jumoTemp1.get()[i];
+        jumoTempData[4*i + 1] = jumoTemp2.get()[i];
+        jumoTempData[4*i + 2] = jumoTemp3.get()[i];
+        jumoTempData[4*i + 3] = jumoTemp4.get()[i];
+    }
+    result.jumoTemp = Aurora::Matrix::New(jumoTempData, 4, jumoTempNum);//输出
+    bool tasTempUsed = false;
+
+    //Tas Temperature
+    if(reconParams::useTASTempComp)
+    {
+        //todo TemperatureModel4D字段不存在
+        //  try
+        // % load temp from files.tempTASComp (standard file, including the preprocessed/corrected temperature data)
+        // % and extract infos, TemperatureModel4D and TASTemperature
+        // TASTemp = loadTASTemperaturesProcessed(path, files.tempTASComp);
+
+        // % extract TemperatureModel4D
+        // if isfield(TASTemp, 'TemperatureModel4D')
+        //     temp.TemperatureModel4D = TASTemp.TemperatureModel4D;
+        // end
+
+        // % extract TASTemperatures
+        // temp.TASTemperature = TASTemp.TASTemperature;
+        // catch
+        //     tasTempCompUsed = 0;
+        // end
+    }
+
+    if (!tasTempUsed)
+    {
+        size_t tasTempLength = aParser->getTemperatureData().getTasTemperature()[0].getLength();
+        float* fromTasTempData = aParser->getTemperatureData().getTasTemperature()[0].get();
+        double* tasTempData = new double[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);
+    }
+
+    std::string speedMethod = "marczak";
+    if (speedMethod == "marczak")
+    {
+        double* kData = new double[6] {2.78786e-9, -1.398845e-6, 3.287156e-4, -5.799136e-2, 5.038813, 1.402385e3};
+        Matrix k = Matrix::fromRawData(kData, 6);
+        result.expectedSOSWater = polyval(k, result.expectedTemp);
+    }
+    else if(speedMethod == "mader")
+    {
+        double* kData = new double[6] {3.14643e-9, -1.478e-6, 0.000334199, -0.0580852, 5.03711, 1402.39};
+        Matrix k = Matrix::fromRawData(kData, 6);
+        result.expectedSOSWater = polyval(k, result.expectedTemp);
+    }
+    else if(speedMethod == "jan")
+    {
+        double speedData = 1557 - 0.0245 * pow((74 - result.expectedTemp.getData()[0]), 2);
+        result.expectedSOSWater = Matrix::fromRawData(new double[1] {speedData}, 1);
+    }
+
+    return result;
+}
+
+//已验证 完全正确
+CEInfo Recon::getCEInfo(Parser* aParser, const MeasurementInfo aInfo)
+{
+    CEInfo result;
+    result.measuredCEUsed = reconParams::useCEMeasured;
+    Matrix ce;
+    auto ceMeasuredPointer = aParser->getAllCEMeasuredData();
+    if (ceMeasuredPointer.isNull())
+    {
+        result.measuredCEAvailable = false;
+        result.measuredCEUsed = false;
+    }
+    else
+    {
+        result.measuredCEAvailable = true;
+        size_t ceMeasuredLength = ceMeasuredPointer.getAscanDataLength();
+        short* ceMeasuredData = ceMeasuredPointer.get();
+        if(aParser->getMetaData().getDataType() == "float16")
+        {
+            ce = convertfp16tofloat(ceMeasuredPointer.get(), ceMeasuredLength, ceMeasuredPointer.getSize());
+        }
+        else
+        {
+            double* ceData = new double[ceMeasuredLength];
+            std::copy(ceMeasuredData, ceMeasuredData + ceMeasuredLength, ceData);
+            ce = Matrix::fromRawData(ceData, ceMeasuredLength, ceMeasuredPointer.getSize());
+        }
+        if(aInfo.Bandpassundersampling)
+        {
+            // todo 2代变量
+            // ceOut(measInfo.NumberSamples, :) = 0; % expand to 1000 could lead to problems in future????
+            // ceOut = reconstructBandpasssubsampling(ceOut, reconstructionFreq, ce.CE_SF);
+        }
+        if (ce.getDimSize(0) < aInfo.numberSamples)
+        {
+            // todo 现有逻辑没有这种可能
+            // indInsert = size(ceOut, 1) + 1:measInfo.NumberSamples;
+            // ceOut(indInsert, :) = repmat(ceOut(size(ceOut, 1), :), length(indInsert), 1); % expand to 1000 could lead to problems in future????
+        }
+        if (reflectParams::removeDCOffset)
+        {
+            ce = ce - repmat(mean(ce),ce.getDimSize(0),1);
+        }
+        result.ce = ce;
+        result.ceOffSet = aParser->getCEData().getCEOffset();
+        result.ce_sf = aInfo.sampleRate;
+        if (aInfo.Hardware == "USCT3dv3")
+        {
+            size_t size = aParser->getMetaData().getTasIndices().getLength();
+            uint8_t* fromData = aParser->getMetaData().getTasIndices().get();
+            double* tasIndicesData = new double[size];
+            std::copy(fromData, fromData + size, tasIndicesData);
+            result.tasIndices = Matrix::fromRawData(tasIndicesData, size);
+
+            size = aParser->getMetaData().getReceiverIndices().getLength();
+            fromData = aParser->getMetaData().getReceiverIndices().get();
+            double* receiverIndicesData = new double[size];
+            std::copy(fromData, fromData + size, receiverIndicesData);
+            result.receiverIndices = Matrix::fromRawData(receiverIndicesData, size);
+        }
+    }
+
+
+    result.ceRefOffSet = aParser->getCEData().getCEOffset();
+    result.ceRef_sf = aParser->getCEData().getCE_SF();
+    auto ceRefFromData = aParser->getCEData().getCE();
+    if (ceRefFromData.isNull())
+    {
+        result.ceAvailable = false;
+    }
+    else
+    {
+        result.ceAvailable = true;
+        size_t ceRefSize = ceRefFromData.getLength();
+        double* ceRefData = new double[ceRefSize];
+        std::copy(ceRefFromData.get(), ceRefFromData.get() + ceRefSize, ceRefData);
+        Matrix ceRefCE = Matrix::fromRawData(ceRefData, ceRefSize);
+        result.ceRef = preprocessCE(ceRefCE, result.ceRef_sf, reflectParams::aScanReconstructionFrequency, aInfo.expectedAScanLength);
+    }
+    result.offsetFilterEnabled = reconParams::offsetFilterEnabled;
+    result.offsetFilterDisabled = reconParams::offsetFilterDisabled;
+
+    return result;
+}

src/common/getMeasurementMetaData.h

@@ -0,0 +1,64 @@
+#ifndef GET_MEASUREMENT_METADATA_H
+#define GET_MEASUREMENT_METADATA_H
+
+#include "Data/TemperatureData.h"
+#include "Matrix.h"
+
+class Parser;
+
+namespace Recon
+{
+    struct MeasurementInfo
+    {
+        unsigned int numberSamples;
+        bool Bandpassundersampling;
+        unsigned int sampleRate;
+        std::string ascanDataType;
+        std::string Hardware;
+        std::string HardwareVersion;
+        int EOffset;
+        int Wavelength;
+        unsigned int expectedAScanLength;
+        std::string rootMeasUniqueID;
+    };
+
+    struct CEInfo
+    {
+        Aurora::Matrix ce;
+        double ceOffSet;
+        double ce_sf;
+        Aurora::Matrix tasIndices;
+        Aurora::Matrix receiverIndices;
+        Aurora::Matrix ceRef;
+        double ceRefOffSet;
+        double ceRef_sf;
+        bool measuredCEUsed;
+        bool measuredCEAvailable;
+        bool ceAvailable;
+        double offsetFilterEnabled;
+        double offsetFilterDisabled;
+    };
+
+    struct TempInfo
+    {
+        Aurora::Matrix jumoTemp;
+        Aurora::Matrix tasTemperature;
+        Aurora::Matrix expectedTemp;
+        Aurora::Matrix expectedSOSWater;
+    };
+
+    struct TransFormInfo
+    {
+        Aurora::Matrix rotationMatrix;
+        Aurora::Matrix motorPos;
+    };
+
+    //无用函数
+    Aurora::Matrix getAvailableMotorPositions(Parser* aParser);
+    MeasurementInfo loadMeasurementInfos(Parser* aParser);
+    TransFormInfo getTransformationMatrix(Parser* aParser, const Aurora::Matrix& aMotorPosList);
+    TempInfo getTemperatureInfo(Parser* aParser, double aNumTas);
+    CEInfo getCEInfo(Parser* aParser, const MeasurementInfo aInfo);
+
+}
+#endif // !GET_MEASUREMENT_METADATA_H

src/config/config.h

@@ -0,0 +1,152 @@
+#ifndef RECON_CONFIG_H
+#define RECON_CONFIG_H
+
+#include <string>
+#include "Matrix.h"
+
+namespace Recon
+{
+    namespace reconParams
+    {
+        //reconParams.measurementInfo.ce
+        static bool useCEMeasured = true;
+        static int removeOutliersFromCEMeasured = 1;
+        static double offsetFilterEnabled = 6.9e-6;
+        static double offsetFilterDisabled = 1.2e-6;
+        //reconParams.measurementInfo.temp
+        static bool useTASTempComp = true;
+        static bool correctTASTemp = 1;
+        static //reconParams.hardwareSelection
+        int gpuSelectionList[8] = {0,1,2,3,4,5,6,7};
+        static //reconParams.dataInfo
+        int expectedAScanDataLength = 4000;
+    }
+
+    namespace reflectParams
+    {
+        //reflectParams.dataSelection
+        static double senderTasList[128];//1~128
+        static double senderElementList[18];//1~18
+        static double receiverTasList[128];//1~128
+        static double receiverElementList[18];//1~18
+        static Aurora::Matrix motorPos = Aurora::Matrix::fromRawData(new double[2]{1,2}, 2);
+        static int constricReflectionAngles = 1;
+        static int angleLowerLimit = 45;
+        static int angleUpperLimit = 360;
+        static int findDefects = 1;
+        static int epsilon = 10;
+        //reflectParams.dataBlocking
+        static int mpSize = 1;
+        static int senderTASSize = 2;
+        static int senderElementSize = 18;
+        //reflectParams.qualityCheck
+        static int qualityCheck = 1;
+        static double warningThreshold = 0.5;
+        static double errorThreshold = 0.1;
+        //reflectParams.dataPreparation
+        static int version = 2;
+        static int aScanReconstructionFrequency = 10000000;
+        static double offsetElectronic =  5.2e-7;
+        static bool removeDCOffset = true;
+        static int expectedAScanDataLength = 4000;
+        //reflectParams.imageInfos
+        static int pixelResolutionX = 400;
+        static int pixelResolutionY = 400;
+        static double imageStartpoint[3] = {-0.18,-0.18,-0.22};
+        static double imageEndpoint[3] = {0.18,0.18,0.02};
+        //reflectParams.signalProcessing
+        static int useOptPulse = 1;
+        static int optPulseFactor = 48;
+        static int expectedPulseLength = 90;
+        static int limitNumPulsesTo = 100;
+        static int normalizePeaks = 1;
+        static int removeTransmissionSignal = 1;
+        static int suppressSameHead = 1;
+        static int suppressSameHeadLength = 1500;
+        static int useCorrelation = 1;
+        static int matchedFilterCeAScan = 1;
+        static int windowLength = 10;
+        static int numThreads = 30;
+        static int expectedUSSpeedRange[2] = {1420,1600};
+        //reflectParams.transmissionCorrection
+        static int soundSpeedCorrection = 1;
+        static int attenuationCorrection = 1;
+        static int saveTransmInReflCoords = 1;
+        static double resolutionTransmMap = 0.005;
+        //reflectParams.saft
+        static int debugMode = 0;
+        static int blockSizeReco = 60000;
+        static int medianWindowSize = 1;
+        static int saftVariant[6] = {1,1,1,1,0,0};
+        static int useAscanIndex = 1;
+        static int attenuationCorrectionLimit = 20;
+        static int blockDimXYZ[3] = {16,16,1};
+        //reflectParams
+        static int VERSION_MATLAB_MAJOR = 9;
+        static int OS_UNIX = 0;
+        static bool runReflectionReco = true;
+    }
+
+    namespace transParams
+    {
+        //transParams.dataSelection
+        static int verbose = 1;
+        static int saveRecon = 1;
+        static int saveDebugInfomation = 1;
+        //transParams.dataSelection
+        static double senderTasList[128];//1~128
+        static double senderElementList[18];//1~18
+        static double receiverTasList[128];//1~128
+        static double receiverElementList[18];//1~18
+        static Aurora::Matrix motorPos = Aurora::Matrix::fromRawData(new double[2]{1,2}, 2);
+        static int filterSensitivity = 1;
+        static double sensFilter = 0.3;
+        static int angleLowerLimit = 120;
+        static int angleUpperLimit = 180;
+        static int applyCalib = 1;
+        static int snrThreshold = 10;
+        static int calibReferenceTas = 1;
+        static int calibReferenceMotorPosition = 1;
+        //transParams.qualityCheck
+        static int qualityCheck = 1;
+        static double warningThreshold = 0.5;
+        static double errorThreshold = 0.1;
+        //transParams.dataBlocking
+        static int mpSize = 1;
+        static int senderTASSize = 2;
+        static int senderElementSize = 18;
+        //transParams.dataPreparation
+        static int numPixelXY = 128;
+        static double offsetElectronic = 5.2e-7;
+        static int aScanReconstructionFrequency = 10000000;
+        static int minTemperature = 15;
+        static int maxTemperature = 40;
+        //transParams.detection
+        static int forceRedetect = 0;
+        static int saveDetection = 0;
+        static int version = 1;
+        static int useTimeWindowing = 1;
+        static int gaussWindow = 0;
+        static int outlierOnTasDetection = 0;
+        static int resampleFactor = 1;
+        static int nThreads = 8;
+        static int minSpeedOfSound = 1450;
+        static int maxSpeedOfSound = 1550;
+        static int detectionWindowSOS = 1;
+        static int detectionWindowATT = 50;
+        //transParams.rayTracing
+        static int bentReconstruction = 0;
+        static int bresenham = 1;
+        static int bentMethod = 1;
+        static int bentTol = 1;
+        static int bentIter = 1;
+        //transParams.solver
+        static std::string name = "TVAL3";
+        static int maxIter = 50;
+        static int muValues = 100;
+        static int betaValues = 1;
+        static bool runTransmissionReco = true;
+    }
+}
+
+#endif //RECON_CONFIG_H