diff --git a/src/reflectionReconstruction/preprocessData/resampleTransmissionVolume.cpp b/src/reflectionReconstruction/preprocessData/resampleTransmissionVolume.cpp new file mode 100644 index 0000000..f8cd2af --- /dev/null +++ b/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 +#include + +using namespace Aurora; +using namespace Recon; + +namespace +{ + Matrix createVectorMatrix(double aStartValue, double aStepValue, double aEndValue) + { + std::vector 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 + +#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