Add buildMatrix(untested)

src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp

@@ -1,11 +1,21 @@
 #include "buildMatrix.h"
+
+#include <cstddef>
+#include <cstdio>
+#include <iostream>
+#include <ostream>
+#include <sys/types.h>
+
 #include "Function1D.h"
 #include "Function2D.h"
 #include "Bresenham.h"
 #include "Function3D.h"
 #include "FMM.h"
-#include <iostream>
-#include <sys/types.h>
+#include "Matrix.h"
+#include "Sparse.h"
+#include "common/common.h"
+#include "config/config.h"
+
 using namespace Aurora;
 namespace Recon
 {
@@ -231,11 +241,109 @@ namespace Recon
         auto gs = precomputeGradients(T);
         //3D
         if(map.getDimSize(2)>1){
-            //TODO:unfinish need function stream3
+            //TODO:unfinish need function stream2
         }
         else{
             //TODO:unfinish need function stream3
         }
         return Matrix();
     }
+
+    BuildMatrixResult buildMatrix(const Aurora::Matrix &senderList,
+                                  const Aurora::Matrix &receiverList,
+                                  Aurora::Matrix& res,
+                                  const Aurora::Matrix &dims, bool BENT,
+                                  const Aurora::Matrix &potentialMap){
+
+        BuildMatrixResult result;
+        size_t numDims = senderList.getDimSize(0);
+        size_t nTotalRays = senderList.getDimSize(1);
+        size_t cnt = 0;
+        size_t coeffIndex = 0;
+        size_t rayCount = 0;
+        size_t safesize = nTotalRays * dims[0] * 2;
+        if (Recon::transParams::saveDebugInfomation){
+            result.hitmap = zeros(dims[0],dims[1],dims[2]);
+        }
+        auto i = zeros(1, safesize);
+        auto j = zeros(1, safesize);
+        auto s = zeros(1, safesize);
+        Matrix path, weighting;
+        int pathLenDisc;
+        bool warnflag = false;
+
+        for (size_t ray = 0; ray < nTotalRays; ray++)
+        {
+            rayCount++;
+            auto startPt = transpose(senderList($,ray).toMatrix());
+            auto endPt = transpose(receiverList($,ray).toMatrix());
+            if (!BENT){
+                if (Recon::transParams::bresenham)
+                {
+                    auto trace = traceStraightRayBresenham(startPt, endPt, res);
+                    path = trace.path;
+                    weighting = trace.weighting;
+                    pathLenDisc = trace.pathLen;
+                }
+                else{
+                    auto trace = traceStraightRay(startPt, endPt, res,dims);
+                    path = trace.path;
+                    weighting = trace.weighting;
+                    pathLenDisc = trace.pathLen;
+                }
+            }
+            else{
+                switch (Recon::transParams::bentMethod) {
+                    case 1:
+                    default:
+                    //TODO:暂时未实现，这里虽然有switch但实际就只有一种。
+                    break;
+                }
+            }
+            Matrix linearIndices =
+                dims.getDimSize(1) == 2
+                    ? convertToLinearIndices(dims, path.block(1, 0, 1))
+                    : convertToLinearIndices(dims, path);
+            if (Recon::transParams::saveDebugInfomation){
+                for (size_t i = 0; i < linearIndices.getDataSize(); i++)
+                {
+                    result.hitmap[linearIndices[i]]+=1;
+                }
+            }
+            printf("Progress: %f (%zu of %zu\r\n)",(double)rayCount*100/(double)nTotalRays,rayCount,nTotalRays);
+            cnt = cnt + pathLenDisc;
+
+            if (cnt < safesize)
+            {
+                i.setBlockValue(1, coeffIndex, cnt, rayCount-1);
+                j.setBlock(1,coeffIndex,cnt,linearIndices);
+                if (weighting.isScalar()){
+                    s.setBlockValue(1, coeffIndex, cnt,weighting.getScalar());
+                }
+                else{
+                    s.setBlock(1, coeffIndex, cnt, weighting);
+                }
+            }
+            else{
+                if(!warnflag)
+                {
+                    warnflag = true;
+                    std::cerr<<"Safesize exceeded, dynamically growing arrays"<<std::endl;
+                }
+                auto rayCountM = ones(1,1);
+                rayCountM[0] = rayCount;
+                i = horzcat(i, repmat(rayCountM,1,pathLenDisc));
+                j = horzcat(i, transpose(linearIndices));
+                if (weighting.isScalar()){
+                    s = horzcat(s, repmat(weighting,1,pathLenDisc));
+                }
+                else{
+                    s = horzcat(s, weighting);
+                }
+            }
+            coeffIndex = cnt;
+        }
+        result.M = Sparse(i.block(0,0,coeffIndex-1), j.block(0,0,coeffIndex-1),s.block(0,0,coeffIndex-1),rayCount,prod(dims).getScalar());
+        return result;
+    }
 }

src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h

@@ -1,6 +1,7 @@
 #ifndef __TRANS__BUILDMATRIX_H__
 #define __TRANS__BUILDMATRIX_H__
 #include "Matrix.h"
+#include "Sparse.h"
 namespace Recon {
     struct TraceStraightRayResult{
         Aurora::Matrix path;
@@ -13,6 +14,11 @@ namespace Recon {
         Aurora::Matrix ds;
     };
 
+    struct BuildMatrixResult{
+        Aurora::Sparse M;
+        Aurora::Matrix hitmap;
+    };
+
     Aurora::Matrix getPixelLengthApproximation(const Aurora::Matrix &startPt,
                                     const Aurora::Matrix &endPt,
                                     const Aurora::Matrix &res, int pathLenDisc);
@@ -31,5 +37,11 @@ namespace Recon {
                                 const Aurora::Matrix &p2,
                                 Aurora::Matrix &discretization);
 
+    BuildMatrixResult buildMatrix(const Aurora::Matrix &senderList,
+                                  const Aurora::Matrix &receiverList,
+                                  Aurora::Matrix& res,
+                                  const Aurora::Matrix &dims, bool BENT,
+                                  const Aurora::Matrix &potentialMap);
+
     }  // namespace Recon
 #endif // __BUILDMATRIX_H__