Add build Matrix function and test to

transmission reconstruction.
2023-05-25 16:19:34 +08:00
parent e98e301ec5
commit 33acaab4fe
5 changed files with 477 additions and 0 deletions
--- a/src/transmissionReconstruction/reconstruction/buildMatrix/FMM.cpp
+++ b/src/transmissionReconstruction/reconstruction/buildMatrix/FMM.cpp
@@ -0,0 +1,100 @@
+#include "FMM.h"
+
+#include <algorithm>
+#include <cstdio>
+#include <iterator>
+#include <sys/types.h>
+#include <tuple>
+#include <vector>
+
+
+#include "Function2D.h"
+#include "Function3D.h"
+#include "Matrix.h"
+#include "Bresenham.h"
+#include "DGradient.h"
+#include "eikonal.h"
+
+
+using namespace Aurora;
+typedef std::tuple<uint,uint,uint> PathRow;
+bool comparePathRow(const PathRow& a,const PathRow& b){
+    if (std::get<0>(a)<std::get<0>(b)) return true;
+    if (std::get<0>(a)==std::get<0>(b))
+    {
+        if(std::get<1>(a)<std::get<1>(b))return true;
+        if (std::get<1>(a)==std::get<1>(b))
+        {
+            return std::get<2>(a)<std::get<2>(b);
+        }
+    }
+    return false;
+}
+namespace Recon
+{
+    void correctPath(Matrix &aMPath, const Matrix &aMStartPt,
+                    const Matrix &aMEndPt) {
+        if (sum(aMStartPt == aMEndPt, Aurora::All).getScalar() ==aMStartPt.getDataSize())return;
+        if (aMPath.isNull() || aMPath.isScalar() || aMPath.isVector()){
+            aMPath = aMEndPt;
+        }
+        int numDims = aMPath.getDimSize(1);
+        auto endPath = aMPath(aMPath.getDimSize(0) - 1, $).toMatrix();
+        size_t size = 0;
+        auto subPathEnd = transpose(Recon::b3dMexDouble(endPath, aMStartPt));
+        std::vector<PathRow> paths;
+        for (size_t i = 0; i < aMPath.getDimSize(0); i++) {
+            uint value3 = numDims > 2 ? aMPath[i + aMPath.getDimSize(0)*2] : 0;
+            paths.push_back({aMPath[i], aMPath[i + aMPath.getDimSize(0)], value3});
+        }
+        for (size_t i = 1; i < subPathEnd.getDimSize(0); i++) {
+            uint value3 =
+                numDims > 2 ? subPathEnd[i + subPathEnd.getDimSize(0) * 2] : 0;
+            paths.push_back(
+                {subPathEnd[i], subPathEnd[i + subPathEnd.getDimSize(0)], value3});
+        }
+        std::sort(paths.begin(), paths.end(), comparePathRow);
+        std::vector<PathRow> path2;
+        std::unique_copy(paths.begin(),paths.end(),std::back_inserter(path2),[](const PathRow& a, const PathRow &b){
+            return std::get<0>(a)==std::get<0>(b) && std::get<1>(a)==std::get<1>(b) && std::get<2>(a)==std::get<2>(b);
+        });
+
+        auto result = zeros(path2.size(),numDims);
+        for (size_t i = 0; i < path2.size(); i++) {
+            printf("path2 index: %d, value: %d, %d, %d\r\n", (int)i,std::get<0>(path2[i]),std::get<1>(path2[i]),std::get<2>(path2[i]));
+            result[i] = std::get<0>(path2[i]);
+            result[i + path2.size()] = std::get<1>(path2[i]);
+            if (numDims > 2)
+                result[i + 2 * path2.size()] = std::get<2>(path2[i]);
+        }
+        aMPath = result;
+    }
+
+    GradientsResult precomputeGradients(Aurora::Matrix& aVa)
+    {
+        GradientsResult result;
+        result.gx = DGradient(aVa, 1, 1);
+        result.gy = DGradient(aVa, 1, 2);
+        result.gz = DGradient(aVa, 1, 3);
+        return result;
+    }
+
+
+    Aurora::Matrix eikonalMex(const Aurora::Matrix& volume, const Aurora::Matrix& point, int gpuSelection){
+        int dims[3]{volume.getDimSize(0), volume.getDimSize(1), volume.getDimSize(2)};
+        // double * result = eikonal(volume.getData(),dims,point.getData(),gpuSelection);
+        double* result=nullptr;
+        return Matrix::fromRawData(result, dims[0], dims[1], dims[2]);
+    }
+
+
+    Aurora::Matrix callFastMarchingMethod(const Aurora::Matrix& map,const Aurora::Matrix& point, const Aurora::Matrix& gpuList)
+    {
+        //TODO:暂时不考虑try catch跳转
+        return eikonalMex(map, point, (int)gpuList[0]);
+    }
+}
+
+
+
+
--- a/src/transmissionReconstruction/reconstruction/buildMatrix/FMM.h
+++ b/src/transmissionReconstruction/reconstruction/buildMatrix/FMM.h
@@ -0,0 +1,23 @@
+#ifndef __RECONSTRUCTION__FMM_H__
+#define __RECONSTRUCTION__FMM_H__
+#include "Matrix.h"
+namespace Recon{
+    struct GradientsResult{
+        Aurora::Matrix gx;
+        Aurora::Matrix gy;
+        Aurora::Matrix gz;
+
+    };
+    void correctPath(Aurora::Matrix& aMPath,const Aurora::Matrix& aMStartPt, const Aurora::Matrix& aMEndPt);
+
+    /**
+     * @brief Computes gradient along dimensions (x, y, z) for given map
+     * 
+     * @param aVa  real double array
+     * @return GradientsResult 
+     */
+    GradientsResult precomputeGradients(Aurora::Matrix& aVa);
+
+    Aurora::Matrix callFastMarchingMethod(const Aurora::Matrix& map, const Aurora::Matrix& point, const Aurora::Matrix& gpuList);
+}
+#endif // __FMM_H__
--- a/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp
+++ b/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp
@@ -0,0 +1,241 @@
+#include "buildMatrix.h"
+#include "Function1D.h"
+#include "Function2D.h"
+#include "Bresenham.h"
+#include "Function3D.h"
+#include "FMM.h"
+#include <iostream>
+#include <sys/types.h>
+using namespace Aurora;
+namespace Recon
+{
+    Aurora::Matrix getPixelLengthApproximation(const Aurora::Matrix& startPt, const Aurora::Matrix& endPt, const Aurora::Matrix& res, int pathLenDisc)
+    {
+        auto dir = endPt - startPt; 
+        auto pathLenReal = sqrt(Aurora::sum((dir * res)^2));
+
+        return pathLenReal / pathLenDisc;
+    }
+    
+    TraceStraightRayResult traceStraightRayBresenham(const Aurora::Matrix &startPt,
+                                  const Aurora::Matrix &endPt,
+                                  const Aurora::Matrix &res)
+    {
+        TraceStraightRayResult result;
+        auto path = transpose(b3dMexDouble(startPt, endPt));
+        result.pathLen = path.getDimSize(0);
+        uint *temp = new uint[path.getDataSize()]{0};
+        //uint32(path)
+        std::copy(path.getData(),path.getData()+path.getDataSize(),temp);
+        std::copy(temp,temp+path.getDataSize(),path.getData()); 
+        delete [] temp;
+        result.path = path;
+        result.weighting = getPixelLengthApproximation(startPt, endPt, res, result.pathLen);
+        return result;
+    }
+
+    int sign(double v){
+        return v==0.0?0:(v<0.0?-1:1);
+    }
+
+    Matrix correctPathToImgDimensions(const Matrix&path, const Matrix&dims)
+    {
+        auto result= path*(path>=1)+(path<1);
+        auto newDims = dims;
+        newDims.forceReshape(1, newDims.getDataSize(), 1);
+        auto repDims = repmat(newDims, result.getDimSize(0), 1);
+        auto invalideValues1 = result <= repDims; 
+        auto invalideValues2 = result > repDims; 
+        result = result*invalideValues1+invalideValues2*repDims;
+        return result;
+    }
+
+    TraceLineResult traceLine3D(const Aurora::Matrix &p1,
+                                    const Aurora::Matrix &p2,
+                                    Aurora::Matrix &discretization)
+    {
+        TraceLineResult result;
+        uint MAX_PATH_LEN = 1500;
+        double EPS = 1e-5;
+        double INFTY = 999999;
+
+        if (discretization.isNull())
+        {
+            discretization = Matrix::fromRawData(new double[3]{1,1,1},1,3);
+        }
+
+        auto dir = p2 - p1;
+        dir = dir/norm(dir,Aurora::Norm2);
+
+        double ta = dir[1] / dir[0];
+
+        double tb = dir[2] / dir[0];
+
+        double tc = dir[1] / dir[2];
+
+        auto pts = zeros( MAX_PATH_LEN, 3);
+        auto ds = zeros( 1, MAX_PATH_LEN);
+        pts( 0, $) = p1;
+        ds[0] = 0;
+        double cnt = 1;
+
+        auto curpt = p1;
+        int sgx = sign( dir[0]);
+        int sgy = sign( dir[1]);
+        int sgz = sign( dir[2]);
+
+        if (std::abs( dir[0]) < EPS)sgx = 0;
+        if (std::abs( dir[1]) < EPS)sgy = 0;
+        if (std::abs( dir[2]) < EPS)sgz = 0;
+
+        while(std::abs( norm( curpt - p2,Norm2)) > EPS){
+            if (cnt > MAX_PATH_LEN){
+                std::cerr<<"cnt > MAX_PATH_LEN="<<MAX_PATH_LEN<<std::endl;
+                return result;
+            }
+            cnt = cnt + 1;
+	        int cc = 0;
+	        if (sgx > 0 ){
+                if(floor( curpt[0] + EPS) == floor( p2[0]))cc++;
+            } 
+	        else if(sgx < 0 ) {
+                if (ceil( curpt[0] - EPS) == ceil( p2[0]))cc++;
+            }else{
+                cc++;
+            } 
+
+            if (sgy > 0){
+                if(floor( curpt[1] + EPS) == floor( p2[1]))cc++;
+            } else if(sgy < 0 ) {
+                if (ceil( curpt[1] - EPS) == ceil( p2[1]))cc++;
+            }else {
+                cc++;
+            }
+            if (sgz > 0 ){
+                if(floor( curpt[2] + EPS) == floor( p2[2]))cc++;
+            } else if(sgz < 0 ) {
+                if (ceil( curpt[2] - EPS) == ceil( p2[2]))cc++;
+            }else {
+                cc++;
+            }
+	
+            if (cc == 3)
+            {
+                pts(cnt-1, $) = p2;
+                ds(cnt-1) = norm( ( curpt - p2) * discretization, Norm2);
+                break;
+            }
+
+            auto a = INFTY;
+            if (sgx > 0)
+            {
+                a = ceil( curpt[0]+ EPS) - curpt[0];
+            }
+            else if( sgx < 0)
+            {
+                a = floor( curpt[0] - EPS) - curpt[0];
+            }
+            double b = INFTY;
+            if (sgy > 0)
+            {
+                b = ceil( curpt[1]+ EPS) - curpt[1];
+            }
+            else if( sgy < 0)
+            {
+                b = floor( curpt[1] - EPS) - curpt[1];
+            }
+            auto c = INFTY;
+            if (sgz > 0)
+            {
+                c = ceil( curpt[2]+ EPS) - curpt[2];
+            }
+            else if( sgz < 0)
+            {
+                c = floor( curpt[2] - EPS) - curpt[2];
+            }
+	
+            double y1 = a * ta;
+            double z1 = a * tb;
+            
+            double x2 = b / ta;
+            double z2 = b / tc;
+
+            double x3 = c / tb;
+            double y3 = c * tc;
+
+	        auto v = transpose(abs(Matrix::fromRawData(new double[9]{a, y1, z1, x2, b, z2, x3, y3, c},3,3)));
+            Matrix idx;
+	        sortrows(v,&idx);
+	        int ix = (int)idx[0];
+            auto nextpt = zeros(1,3);
+            if (ix == 0)
+            {
+                nextpt[0] = curpt[0] + a;
+                nextpt[1] = curpt[1] + y1;
+                nextpt[2] = curpt[2] + z1;
+            }
+            else if (ix == 1){
+                nextpt[0] = curpt[0] + x2;
+                nextpt[1] = curpt[1] + b;
+                nextpt[2] = curpt[2] + z2;
+            }
+            else{
+                nextpt[0] = curpt[0] + x3;
+                nextpt[1] = curpt[1] + y3;
+                nextpt[2] = curpt[2] + c;
+            }
+            pts( cnt-1, $) = nextpt;
+            ds[cnt-1] = norm( ( curpt - nextpt) * discretization, Norm2);
+            curpt = nextpt;
+        }
+        result.path = pts.block(0, 0, cnt-1);
+        result.ds = ds.block(1, 0, cnt-1);
+        return result;
+    }
+    
+    TraceStraightRayResult traceStraightRay(const Aurora::Matrix &startPt,
+                                                const Aurora::Matrix &endPt,
+                                                Aurora::Matrix &res,
+                                                const Aurora::Matrix & dims)
+    {
+        TraceStraightRayResult result;
+        int numDims = Aurora::size(startPt, 2);
+        if (numDims==2){
+
+        }
+        else if (numDims!=3){
+            std::cerr<<"Number of values "<<numDims<<" not supported. Needs to be coordinates for 2D or 3D"<<std::endl;
+            return result;
+        }
+        auto traceLineResult = traceLine3D(startPt, endPt, res);
+
+        auto path = traceLineResult.path.block(1, 0,numDims-1);
+
+        uint *temp = new uint[path.getDataSize()]{0};
+        //uint32(path)
+        std::copy(path.getData(),path.getData()+path.getDataSize(),temp);
+        std::copy(temp,temp+path.getDataSize(),path.getData()); 
+        delete [] temp;
+        result.weighting = traceLineResult.ds;
+        result.path = correctPathToImgDimensions(path, dims);
+        result.pathLen = result.path.getDimSize(0);
+        return result;
+    }
+
+    //TODO:unfinished!
+    Aurora::Matrix traceBentRayFM(const Aurora::Matrix &startPoint,
+                                  const Aurora::Matrix &endPoints,
+                                  const Aurora::Matrix &map,
+                                  const Aurora::Matrix &gputList) {
+        auto T = callFastMarchingMethod(map,startPoint,gputList);
+        auto gs = precomputeGradients(T);
+        //3D
+        if(map.getDimSize(2)>1){
+            //TODO:unfinish need function stream3
+        }
+        else{
+            //TODO:unfinish need function stream3
+        }
+        return Matrix();
+    }
+}
--- a/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h
+++ b/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h
@@ -0,0 +1,35 @@
+#ifndef __TRANS__BUILDMATRIX_H__
+#define __TRANS__BUILDMATRIX_H__
+#include "Matrix.h"
+namespace Recon {
+    struct TraceStraightRayResult{
+        Aurora::Matrix path;
+        Aurora::Matrix weighting;
+        int pathLen;
+    };
+    
+    struct TraceLineResult{
+        Aurora::Matrix path;
+        Aurora::Matrix ds;
+    };
+
+    Aurora::Matrix getPixelLengthApproximation(const Aurora::Matrix &startPt,
+                                    const Aurora::Matrix &endPt,
+                                    const Aurora::Matrix &res, int pathLenDisc);
+
+    TraceStraightRayResult
+    traceStraightRayBresenham(const Aurora::Matrix &startPt,
+                              const Aurora::Matrix &endPt,
+                              const Aurora::Matrix &res);
+
+    TraceStraightRayResult traceStraightRay(const Aurora::Matrix &startPt,
+                                            const Aurora::Matrix &endPt,
+                                            Aurora::Matrix &res,
+                                            const Aurora::Matrix & dims);
+
+    TraceLineResult traceLine3D(const Aurora::Matrix &p1,
+                                const Aurora::Matrix &p2,
+                                Aurora::Matrix &discretization);
+
+    }  // namespace Recon
+#endif // __BUILDMATRIX_H__