From 33acaab4fe939ad4283e25443dea83264a19c1b8 Mon Sep 17 00:00:00 2001 From: kradchen Date: Thu, 25 May 2023 16:19:34 +0800 Subject: [PATCH] Add build Matrix function and test to transmission reconstruction. --- .../reconstruction/buildMatrix/FMM.cpp | 100 ++++++++ .../reconstruction/buildMatrix/FMM.h | 23 ++ .../buildMatrix/buildMatrix.cpp | 241 ++++++++++++++++++ .../reconstruction/buildMatrix/buildMatrix.h | 35 +++ test/Reconstruction_Test.cpp | 78 ++++++ 5 files changed, 477 insertions(+) create mode 100644 src/transmissionReconstruction/reconstruction/buildMatrix/FMM.cpp create mode 100644 src/transmissionReconstruction/reconstruction/buildMatrix/FMM.h create mode 100644 src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp create mode 100644 src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h diff --git a/src/transmissionReconstruction/reconstruction/buildMatrix/FMM.cpp b/src/transmissionReconstruction/reconstruction/buildMatrix/FMM.cpp new file mode 100644 index 0000000..63be901 --- /dev/null +++ b/src/transmissionReconstruction/reconstruction/buildMatrix/FMM.cpp @@ -0,0 +1,100 @@ +#include "FMM.h" + +#include +#include +#include +#include +#include +#include + + +#include "Function2D.h" +#include "Function3D.h" +#include "Matrix.h" +#include "Bresenham.h" +#include "DGradient.h" +#include "eikonal.h" + + +using namespace Aurora; +typedef std::tuple PathRow; +bool comparePathRow(const PathRow& a,const PathRow& b){ + if (std::get<0>(a)(b)) return true; + if (std::get<0>(a)==std::get<0>(b)) + { + if(std::get<1>(a)(b))return true; + if (std::get<1>(a)==std::get<1>(b)) + { + return std::get<2>(a)(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 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 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]); + } +} + + + + diff --git a/src/transmissionReconstruction/reconstruction/buildMatrix/FMM.h b/src/transmissionReconstruction/reconstruction/buildMatrix/FMM.h new file mode 100644 index 0000000..ba120b1 --- /dev/null +++ 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__ \ No newline at end of file diff --git a/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp b/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp new file mode 100644 index 0000000..2fc7f8c --- /dev/null +++ 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 +#include +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="< 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 "<1){ + //TODO:unfinish need function stream3 + } + else{ + //TODO:unfinish need function stream3 + } + return Matrix(); + } +} \ No newline at end of file diff --git a/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h b/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h new file mode 100644 index 0000000..1f855be --- /dev/null +++ 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__ \ No newline at end of file diff --git a/test/Reconstruction_Test.cpp b/test/Reconstruction_Test.cpp index 9fce757..a5b85ef 100644 --- a/test/Reconstruction_Test.cpp +++ b/test/Reconstruction_Test.cpp @@ -5,6 +5,8 @@ #include "MatlabReader.h" #include "Matrix.h" #include "transmissionReconstruction/reconstruction/buildMatrix/DGradient.h" +#include "transmissionReconstruction/reconstruction/buildMatrix/FMM.h" +#include "transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h" #include "transmissionReconstruction/reconstruction/reconstruction.h" @@ -85,3 +87,79 @@ TEST_F(Reconstruction_Test, DGradient) { } } +TEST_F(Reconstruction_Test, correctPath) { + auto path = Aurora::Matrix::fromRawData(new double[6]{1,1,10,9,1,1},2,3); + auto startPt = Aurora::Matrix::fromRawData(new double[3]{1, .5, 1}, 1, 3); + auto endPt = Aurora::Matrix::fromRawData(new double[3]{1, 10, 1}, 1, 3); + + Recon::correctPath(path,startPt,endPt); + for (size_t i = 0; i < 10; i++) + { + EXPECT_DOUBLE_EQ(path[i],1.0); + EXPECT_DOUBLE_EQ(path[i+10],1.0+i); + EXPECT_DOUBLE_EQ(path[i+20],1.0); + } +} + +TEST_F(Reconstruction_Test,getPixelLengthApproximation){ + auto startPt = Aurora::Matrix::fromRawData(new double[3]{1, 1, 1}, 1, 3); + auto endPt = Aurora::Matrix::fromRawData(new double[3]{4, 5, 6}, 1, 3); + auto res = Aurora::Matrix::fromRawData(new double[3]{.1, .1, .1}, 1, 3); + auto weight = Recon::getPixelLengthApproximation(startPt, endPt, res, 10); + EXPECT_DOUBLE_AE(weight[0],.0707); +} + +TEST_F(Reconstruction_Test,traceLine3D){ + auto p1 = Aurora::Matrix::fromRawData(new double[3]{1, 10, 12}, 1, 3); + auto p2 = Aurora::Matrix::fromRawData(new double[3]{20, 2, 13}, 1, 3); + auto discretization = Aurora::Matrix::fromRawData(new double[3]{1, 1, 1}, 1, 3); + auto result = Recon::traceLine3D(p1, p2, discretization); + MatlabReader m("/home/krad/TestData/traceLine3D.mat"); + + auto path = m.read("path"); + auto ds = m.read("ds"); + result.path.printf(); + for (size_t i = 0; i < result.path.getDataSize(); i++) + { + EXPECT_DOUBLE_AE(result.path.getData()[i],path[i])<<"index:"<