Add discretizePositions to Transmission

recontruction.

src/transmissionReconstruction/reconstruction/reconstruction.cpp

@@ -2,7 +2,10 @@
 #include "Function.h"
 #include "Function1D.h"
 #include "Function2D.h"
+#include "Function3D.h"
 #include "Matrix.h"
+#include <cmath>
+#include <iostream>
 using namespace Aurora;
 namespace Recon {
     Aurora::Matrix calculateMinimalMaximalTransducerPositions(
@@ -32,4 +35,60 @@ namespace Recon {
         ddims.forceReshape(1, ddims.getDataSize(), 1);
         return  ddims;
     }
+    
+    Aurora::Matrix calculateResolution(const Aurora::Matrix &aVDdims, const Aurora::Matrix &aVDims)
+    {
+        auto numDim = aVDims.getDataSize();
+        auto res = (aVDdims.block(1,numDim,aVDdims.getDataSize()-1) - aVDdims.block(1,0,numDim-1))/ (aVDims - 1);
+
+        if(numDim == 3 && aVDims[2] == 1){
+            res[2] = res[0];
+        }
+        return res;
+    }
+    
+    Aurora::Matrix getDimensions(double aNumPixelXY, const Aurora::Matrix& aVDdims)
+    {
+        int numDim = aVDdims.getDataSize()/2;
+
+        if (!(numDim == 2 || numDim == 3)){
+           std::cerr<<"Inputs does not match requirements."<<std::endl;
+           return Aurora::Matrix();
+        }
+
+        auto dims = ones(1,numDim);
+
+        dims[0] = aNumPixelXY;
+        dims[1] = aNumPixelXY; 
+
+        if (numDim == 3)
+        {
+            dims[2] = std::ceil((aVDdims[5]-aVDdims[2])/((aVDdims[3]-aVDdims[0])/dims[0]));
+            //BUG:发现源代码可能存在无限值，加了std::isfinite验证
+            if(dims[2] == 0 || !std::isfinite(dims[2]))dims[2] = 1;
+        }
+        return dims;
+    }
+    
+    void slownessToSOS(Aurora::Matrix & aVF1, double aSOS_IN_WATER)
+    {
+        aVF1 = 1 / ((aVF1 + 1) / aSOS_IN_WATER);
+    }
+    
+    DiscretizePositionValues discretizePositions(Aurora::Matrix &aVSenderCoordList, Aurora::Matrix &aVReceiverCoordList, double aNumPixelXY)
+    {
+        DiscretizePositionValues result;
+        result.ddims = calculateMinimalMaximalTransducerPositions(aVSenderCoordList, aVReceiverCoordList);
+        result.dims = getDimensions(aNumPixelXY, result.ddims);
+        result.res = calculateResolution(result.ddims, result.dims);
+        int numDim = aVSenderCoordList.getDimSize(0);
+        auto tempDdims = result.ddims.block(0,0,numDim-1);
+        auto tempTransRes = repmat(transpose(result.res),1,2);
+        auto tempDivRes = repmat(transpose(tempDdims / result.res),1,2);
+        aVSenderCoordList = (aVSenderCoordList / tempTransRes) - tempDivRes + 1;
+        aVReceiverCoordList = (aVReceiverCoordList / tempTransRes) - tempDivRes + 1;
+        result.receiverCoordList = aVReceiverCoordList;
+        result.senderCoordList = aVSenderCoordList;
+        return result;
+    }
 } // namespace Recon

src/transmissionReconstruction/reconstruction/reconstruction.h

@@ -2,7 +2,25 @@
 #define _TRANS_RECONSTRUCTION_H__
 #include "Matrix.h"
 namespace Recon {
+struct DiscretizePositionValues{
+   Aurora::Matrix senderCoordList;
+   Aurora::Matrix receiverCoordList;
+   Aurora::Matrix ddims;
+   Aurora::Matrix dims;
+   Aurora::Matrix res;
+};
 Aurora::Matrix calculateMinimalMaximalTransducerPositions(
-    const Aurora::Matrix &aMSenderList, const Aurora::Matrix &aMReceiverList);
+   const Aurora::Matrix &aMSenderList, const Aurora::Matrix &aMReceiverList);
+
+Aurora::Matrix calculateResolution(const Aurora::Matrix &aVDdims, const Aurora::Matrix &aVDims);
+
+Aurora::Matrix getDimensions(double aNumPixelXY, const Aurora::Matrix& ddims);
+
+void slownessToSOS(Aurora::Matrix & aVF1, double aSOS_IN_WATER);
+
+DiscretizePositionValues discretizePositions(Aurora::Matrix &aVSenderCoordList, Aurora::Matrix &aVReceiverCoordList, double aNumPixelXY);
+
+//TODO: exitBent, reconstructArt, saveTransmissionReconstruction.m
+
 }
 #endif // __RECONSTRUCTION_H__

test/Reconstruction_Test.cpp

@@ -41,7 +41,33 @@ TEST_F(Reconstruction_Test, calculateMinimalMaximalTransducerPositions) {
     EXPECT_DOUBLE_EQ(1.0,result.getData()[2]);
     EXPECT_DOUBLE_EQ(1.0,result.getData()[3]);
     EXPECT_DOUBLE_EQ(8,result.getData()[4]);
-    EXPECT_DOUBLE_EQ(4,result.getData()[5]);
+    EXPECT_DOUBLE_EQ(4,result.getData()[5]); 
-
-    
 }
+
+TEST_F(Reconstruction_Test, calculateResolution) {
+    auto ddims   = Aurora::Matrix::fromRawData(new double[6]{-0.1296,-0.1296,0.0185,0.1296,0.1295,0.1682}, 1, 6);
+    auto dims  = Aurora::Matrix::fromRawData(new double[3]{128,128,74}, 1, 3);
+    auto result  = Recon::calculateResolution(ddims, dims);
+    EXPECT_DOUBLE_AE(0.0020,result[0]);
+    EXPECT_DOUBLE_AE(0.0020,result[1]);
+    EXPECT_DOUBLE_AE(0.0021,result[2]);
+}
+
+TEST_F(Reconstruction_Test, getDimensions) {
+    auto ddims = Aurora::Matrix::fromRawData(new double[6]{-0.1296,-0.1296,0.0185,0.1296,0.1295,0.1682}, 1, 6);
+    double numPixelXY = 128;
+    auto result  = Recon::getDimensions(numPixelXY,ddims);
+    EXPECT_DOUBLE_AE(128,result[0]);
+    EXPECT_DOUBLE_AE(128,result[1]);
+    EXPECT_DOUBLE_AE(74,result[2]);
+}
+
+TEST_F(Reconstruction_Test, discretizePositions) {
+    auto senderList  = Aurora::Matrix::fromRawData(new double[6]{1, 2, 3, 1, 2, 4}, 3, 2);
+    auto receiverList  = Aurora::Matrix::fromRawData(new double[6]{1, 8, 3, 1, 2, 1}, 3, 2);
+    double numPixelXY = 128;
+    auto result  = Recon::discretizePositions(senderList,receiverList,numPixelXY);
+    EXPECT_DOUBLE_AE(128,result.receiverCoordList[1]);
+    EXPECT_DOUBLE_AE(1,result.senderCoordList[1]);
+}
+