diff --git a/CMakeLists.txt b/CMakeLists.txt index b247743..a1e4dd2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,7 +22,7 @@ target_link_libraries(UR PUBLIC ${Parser_Libraries}) target_link_libraries(UR PUBLIC URDepends::TransDetection) target_link_libraries(UR PUBLIC URDepends::eikonal) - +# target_link_libraries(UR PUBLIC URDepends::TVALGPU) find_package(GTest REQUIRED) INCLUDE_DIRECTORIES(${GTEST_INCLUDE_DIRS}) @@ -41,4 +41,5 @@ target_link_libraries(UR_Test PUBLIC matio) target_link_libraries(UR_Test PUBLIC ${Parser_Libraries}) target_link_libraries(UR_Test PUBLIC URDepends::TransDetection) target_link_libraries(UR_Test PUBLIC URDepends::eikonal) +# target_link_libraries(UR_Test PUBLIC URDepends::TVALGPU) gtest_discover_tests(UR_Test) \ No newline at end of file diff --git a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp new file mode 100644 index 0000000..b7cebc4 --- /dev/null +++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp @@ -0,0 +1,30 @@ +#include "TVAL.h" + +#include "Matrix.h" +#include "tval3gpu3d.h" +#include +#include + +namespace Recon +{ + Aurora::Matrix callTval3(Aurora::Sparse& M, const Aurora::Matrix& b,const Aurora::Matrix& dims,int device, const TVALOptions& opt) + { + int * xIdxs = new int[M.getColVector().getDataSize()]; + std::copy(M.getColVector().getData(),M.getColVector().getData()+M.getColVector().getDataSize(),xIdxs); + int * yIdxs = new int[M.getRowVector().getDataSize()]; + std::copy(M.getRowVector().getData(),M.getRowVector().getData()+M.getRowVector().getDataSize(),yIdxs); + Aurora::Matrix values = M.getValVector(); + int nz = std::max(M.getColVector().getDataSize(),std::max(M.getRowVector().getDataSize(),M.getValVector().getDataSize())); + float* bData = new float[b.getDataSize()]; + std::copy(b.getData(),b.getData()+b.getDataSize(),bData); + size_t bDims[3]={(size_t)b.getDimSize(0),(size_t)b.getDimSize(1),(size_t)b.getDimSize(2)}; + size_t rdims[3] = {(size_t)dims[0], (size_t)dims[1], (size_t)dims[2]}; + bool pagelocked = false; + // auto result = TVALGPU(xIdxs, yIdxs, values.getData(), M.getM(), M.getN(), nz, bData, bDims, rdims, opt, device, false); + delete [] xIdxs; + delete [] yIdxs; + delete [] bData; + // return Aurora::Matrix::fromRawData(result.data, result.dims[0],result.dims[1],result.dims[2]); + return Aurora::Matrix(); + } +} \ No newline at end of file diff --git a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h new file mode 100644 index 0000000..4800006 --- /dev/null +++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h @@ -0,0 +1,11 @@ +#ifndef __TVAL_H__ +#define __TVAL_H__ +#include "Matrix.h" +#include "Sparse.h" +#include "tvalstruct.h" + +namespace Recon { + Aurora::Matrix callTval3(Aurora::Sparse& M, const Aurora::Matrix& b,const Aurora::Matrix& dims, int device, const struct TVALOptions& options); +} + +#endif // __TVAL_H__ \ No newline at end of file diff --git a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.cpp b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.cpp new file mode 100644 index 0000000..117caa3 --- /dev/null +++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.cpp @@ -0,0 +1,98 @@ +#include "solve.h" + +#include + +#include "Function3D.h" +#include "Matrix.h" + +#include "config/config.h" + +#include "transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h" +#include "tvalstruct.h" + +namespace Recon +{ + struct TVAL3SolverOptions{ + Aurora::Matrix gpuSelectionList; + double TVAL3MU; + double TVAL3MU0; + double TVAL3Beta; + double TVAL3Beta0; + bool nonNeg = false; + }; + + Aurora::Matrix solve( Aurora::Sparse& M, const Aurora::Matrix& b, const Aurora::Matrix& dims, int niter, TVAL3SolverOptions solverOptions){ + if (Recon::transParams::name.empty()){ + Recon::transParams::name = "TVAL3"; + } + if (Recon::transParams::name == "TVAL3") + { + //callTval3 + TVALOptions opt; + opt.bent = false; + opt.tol = 1E-10; + opt.maxit = niter; + opt.TVnorm = 2; + opt.disp = false; + opt.mu0 = solverOptions.TVAL3MU0; + opt.mu = solverOptions.TVAL3MU; + opt.beta = solverOptions.TVAL3Beta; + opt.beta0 = solverOptions.TVAL3Beta0; + int device = (int)solverOptions.gpuSelectionList[0]; + return callTval3(M, b, dims, device, opt); + } + //SART + else{ + //TODO:待实现，先实现默认的TVAL3 + return Aurora::Matrix(); + } + } + + std::vector> solveParameterIterator(Aurora::Sparse M, const Aurora::Matrix &b, + const Aurora::Matrix &dims, bool oneIter, bool nonNeg) + { + if (Recon::transParams::name == "TVAL3"){ + std::vector> result(Recon::transParams::muValues.getDataSize()); + if (Recon::transParams::muValues.isNull()){ + Recon::transParams::muValues = Aurora::ones(1,1); + Recon::transParams::muValues[0] = 24; + } + if (Recon::transParams::betaValues.isNull()){ + Recon::transParams::betaValues = Aurora::ones(1,1); + Recon::transParams::betaValues[0] = 1; + } + if (oneIter){ + auto temp = Aurora::ones(1,1); + temp[0] = Recon::transParams::muValues[0]; + Recon::transParams::muValues = temp; + auto temp2 = Aurora::ones(1,1); + temp2[0] = Recon::transParams::betaValues[0]; + Recon::transParams::betaValues = temp2; + } + TVAL3SolverOptions options; + options.gpuSelectionList = Recon::transParams::gpuSelectionList; + for (size_t i = 0; i < Recon::transParams::muValues.getDataSize(); i++) + { + options.TVAL3MU = Recon::transParams::muValues[i]; + options.TVAL3MU0 = Recon::transParams::muValues[i]; + std::vector solveResult(Recon::transParams::betaValues.getDataSize()); + for (size_t j = 0; j < Recon::transParams::betaValues.getDataSize(); j++) + { + options.TVAL3Beta = Recon::transParams::betaValues[i]; + options.TVAL3Beta0 = Recon::transParams::betaValues[i]; + options.nonNeg = nonNeg; + solveResult[j] = solve(M, b, dims, transParams::maxIter, options); + } + result[i] = solveResult; + } + return result; + } + //SART + else{ + std::vector> result; + //TODO:暂时未实现 + return result; + } + + } +} \ No newline at end of file diff --git a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.h b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.h new file mode 100644 index 0000000..2b097a6 --- /dev/null +++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.h @@ -0,0 +1,12 @@ +#ifndef __SOLVE_H__ +#define __SOLVE_H__ +#include "Matrix.h" +#include "Sparse.h" +#include +namespace Recon { + std::vector> + solveParameterIterator(Aurora::Sparse M, const Aurora::Matrix &b, + const Aurora::Matrix &dims, bool oneIter = true, bool nonNeg = false); +} // namespace Recon + +#endif // __SOLVE_H__ \ No newline at end of file