Add resampleTransmissionVolume and unittest.

src/reflectionReconstruction/preprocessData/resampleTransmissionVolume.cpp

@@ -0,0 +1,111 @@
+#include "resampleTransmissionVolume.h"
+#include "imageExtrapolation.h"
+#include "config/config.h"
+#include "src/common/meshgrid.h"
+
+#include "Matrix.h"
+#include "Function2D.h"
+#include "Function1D.h"
+#include "Function3D.h"
+#include <cstddef>
+#include <vector>
+
+using namespace Aurora;
+using namespace Recon;
+
+namespace 
+{
+    Matrix createVectorMatrix(double aStartValue, double aStepValue, double aEndValue)
+    {
+        std::vector<double> matrixData;
+        double tempValue = aStartValue;
+        matrixData.push_back(tempValue);
+        long long compare1 = round(aEndValue * 10e13);
+        long long compare2 = round(tempValue * 10e13);
+        while(round(tempValue* 10e13) <= compare1)
+        {
+            tempValue += aStepValue;
+            matrixData.push_back(tempValue);
+            compare2 = round(tempValue * 10e14);
+        }
+        matrixData.pop_back();
+
+        return Matrix::copyFromRawData(matrixData.data(), 1, matrixData.size());
+    }
+}
+
+TransmissionVolme Recon::resampleTransmissionVolume(const Aurora::Matrix& aTransVol, const Aurora::Matrix& aMaxPos, 
+                                                 const Aurora::Matrix& aMinPos,const GeometryInfo& aGeom)
+{
+    TransmissionVolme result;
+    double resSOSX = (aMaxPos[0] - aMinPos[0]) / (aTransVol.getDimSize(0) - 1);
+    double resSOSY = (aMaxPos[1] - aMinPos[1]) / (aTransVol.getDimSize(1) - 1);
+    double resSOSZ = (aMaxPos[2] - aMinPos[2]) / (aTransVol.getDimSize(2) - 1);
+
+    Matrix xSOS = createVectorMatrix(aMinPos[0] + 0.5 * resSOSX, resSOSX, aMaxPos[0] + 0.5 * resSOSX);
+    Matrix ySOS = createVectorMatrix(aMinPos[1] + 0.5 * resSOSY, resSOSY, aMaxPos[1] + 0.5 * resSOSY);
+    Matrix zSOS = createVectorMatrix(aMinPos[2] + 0.5 * resSOSZ, resSOSZ, aMaxPos[2] + 0.5 * resSOSZ);
+    xSOS = transpose(xSOS);
+    ySOS = transpose(ySOS);
+    zSOS = transpose(zSOS);
+
+    Matrix minTransdPos = min(aGeom.minEmitter, aGeom.minReceiver);
+    Matrix maxTransdPos = max(aGeom.maxEmitter, aGeom.maxReceiver);
+    double deltaTransMap = (aMaxPos[0] - aMinPos[0]) / (aTransVol.getDimSize(0) - 1);
+    Matrix beginTransMap = transpose(min(transpose(reflectParams::imageStartpoint), minTransdPos)) - 3 * reflectParams::resolutionTransmMap;
+    Matrix endSpeedMap = max(transpose(reflectParams::imageEndpoint), maxTransdPos) + 3 * reflectParams::resolutionTransmMap;
+    endSpeedMap[2] = endSpeedMap[2] + 0.025;    
+
+    Matrix xNew = createVectorMatrix(beginTransMap[0], deltaTransMap, endSpeedMap[0]) + 0.5 * reflectParams::resolutionTransmMap;
+    Matrix yNew = createVectorMatrix(beginTransMap[1], deltaTransMap, endSpeedMap[1]) + 0.5 * reflectParams::resolutionTransmMap;
+    Matrix zNew = createVectorMatrix(beginTransMap[2], deltaTransMap, endSpeedMap[2]) + 0.5 * reflectParams::resolutionTransmMap;
+    xNew = transpose(xNew);
+    yNew = transpose(yNew);
+    zNew = transpose(zNew);
+    
+    Matrix smHelper = meshgridInterp3(ySOS,xSOS,zSOS,aTransVol,yNew,xNew,zNew, InterpnMethod::Spline, 0);
+    //result.transMap = smHelper;
+    smHelper = imageExtrapolation(smHelper, 3);
+    xSOS = xNew;
+    ySOS = yNew;
+    zSOS = zNew;
+    Matrix transVol = smHelper;
+    for(size_t i=0; i<transVol.getDataSize(); ++i)
+    {
+        if(transVol[i] == 0)
+        {
+            transVol[i] = -1e7;
+        }
+    }
+
+    if(reflectParams::resolutionTransmMap != 0)
+    {
+        deltaTransMap = reflectParams::resolutionTransmMap;
+        xNew = createVectorMatrix(beginTransMap[0], deltaTransMap, endSpeedMap[0]) + 0.5 * deltaTransMap;
+        yNew = createVectorMatrix(beginTransMap[1], deltaTransMap, endSpeedMap[1]) + 0.5 * deltaTransMap;
+        zNew = createVectorMatrix(beginTransMap[2], deltaTransMap, endSpeedMap[2]) + 0.5 * deltaTransMap;
+        xNew = transpose(xNew);
+        yNew = transpose(yNew);
+        zNew = transpose(zNew);
+        smHelper = meshgridInterp3(ySOS,xSOS,zSOS,transVol,yNew,xNew,zNew, InterpnMethod::Linear, 0);
+    }
+    else
+    {
+        smHelper = transVol;
+    }
+    
+    for(size_t i=0; i<smHelper.getDataSize(); ++i)
+    {
+        if(smHelper[i] < 0)
+        {
+            smHelper[i] = 0;
+        }
+    }
+
+    smHelper = imageExtrapolation(smHelper, 3);
+    result.transMap = smHelper;
+    result.deltaTransMap = deltaTransMap;
+    result.beginTransMap = beginTransMap;
+
+    return result;
+}                                                 

src/reflectionReconstruction/preprocessData/resampleTransmissionVolume.h

@@ -0,0 +1,21 @@
+#ifndef RESAMPLETRANSMISSIONVOLUME_H
+#define RESAMPLETRANSMISSIONVOLUME_H
+
+#include "Matrix.h"
+#include "src/common/getGeometryInfo.h"
+
+namespace Recon
+{
+    struct TransmissionVolme
+    {
+        Aurora::Matrix transMap;
+        double deltaTransMap;
+        Aurora::Matrix beginTransMap;
+    };
+
+    TransmissionVolme resampleTransmissionVolume(const Aurora::Matrix& aTransVol, const Aurora::Matrix& aMaxPos, 
+                                                 const Aurora::Matrix& aMinPos,const GeometryInfo& aGeom);
+}
+
+
+#endif // RESAMPLETRANSMISSIONVOLUME_H

test/ProcessData_Test.cpp

@@ -0,0 +1,86 @@
+#include <gtest/gtest.h>
+
+#include "Function2D.h"
+#include "Matrix.h"
+
+#include "src/reflectionReconstruction/preprocessData/resampleTransmissionVolume.h"
+#include "src/reflectionReconstruction/preprocessData/precalcImageParameters.h"
+#include "src/reflectionReconstruction/preprocessData/imageExtrapolation.h"
+#include "src/common/getGeometryInfo.h"
+#include "src/config/config.h"
+
+#include "MatlabReader.h"
+
+using namespace Recon;
+
+inline double fourDecimalRound(double src){
+    return round(src*10000.0)/10000.0;
+}
+
+#define EXPECT_DOUBLE_AE(valueA,valueB)\
+    EXPECT_DOUBLE_EQ(fourDecimalRound(valueA),fourDecimalRound(valueB))
+
+class ProcessData_Test : public ::testing::Test {
+protected:
+    static void SetUpCommonTester() {
+
+    }
+
+    static void TearDownTestCase() {
+    }
+
+    void SetUp() {
+    }
+
+    void TearDown() {
+    }
+};
+
+TEST_F(ProcessData_Test, imageExtrapolation) {
+    MatlabReader m("/home/sun/testData/imageExtrapolation.mat");
+
+    auto smHelperInput = m.read("SMHelperInput");
+    auto smHelper = m.read("SMHelper");
+    auto result = imageExtrapolation(smHelperInput, 3);
+    for(int i=0;i<result.getDataSize();++i)
+    {
+        EXPECT_DOUBLE_AE(smHelper[i],result[i]) << i;
+    }
+    
+}
+
+TEST_F(ProcessData_Test, resampleTransmissionVolume) {
+    MatlabReader m("/home/sun/testData/resampleTransmissionVolume.mat");
+
+    auto maxEmitter = m.read("maxEmitter");
+    auto maxPos = m.read("maxPos");
+    auto maxReceiver = m.read("maxReceiver");
+    auto minEmitter = m.read("minEmitter");
+    auto minPos = m.read("minPos");
+    auto minReceiver = m.read("minReceiver");
+    auto transVol = m.read("transVol");
+    auto speedMap3d = m.read("SpeedMap3D");
+    auto beginTransMap = m.read("begin_TransMap");
+    auto deltaTransMap = m.read("delta_TransMap");
+    Recon::GeometryInfo geom;
+    geom.maxEmitter = maxEmitter;
+    geom.minEmitter = minEmitter;
+    geom.maxReceiver = maxReceiver;
+    geom.minReceiver = minReceiver;
+    geom.minSize = m.read("minSize");
+    geom.maxSize = m.read("maxSize");
+    precalcImageParameters(geom);
+    auto r = resampleTransmissionVolume(transVol, maxPos, minPos, geom);
+    for(int i=0;i<r.transMap.getDataSize();++i)
+    {
+        EXPECT_DOUBLE_AE(speedMap3d[i],r.transMap[i]);
+    }
+    for(int i=0;i<r.beginTransMap.getDataSize();++i)
+    {
+        EXPECT_DOUBLE_AE(beginTransMap[i],r.beginTransMap[i]);
+    }
+
+    EXPECT_DOUBLE_AE(deltaTransMap[0],r.deltaTransMap);
+
+}
+