Add solve and TVAL(un finished).

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)

src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.cpp

@@ -0,0 +1,30 @@
+#include "TVAL.h"
+
+#include "Matrix.h"
+#include "tval3gpu3d.h"
+#include <algorithm>
+#include <cstddef>
+
+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();
+    }
+}

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__

src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.cpp

@@ -0,0 +1,98 @@
+#include "solve.h"
+
+#include <cstddef>
+
+#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<std::vector<Aurora::Matrix>> solveParameterIterator(Aurora::Sparse M, const Aurora::Matrix &b,
+                               const Aurora::Matrix &dims, bool oneIter, bool nonNeg)
+    {
+        if (Recon::transParams::name == "TVAL3"){
+            std::vector<std::vector<Aurora::Matrix>> 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<Aurora::Matrix> 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<std::vector<Aurora::Matrix>> result;
+            //TODO:暂时未实现
+            return result;
+        }
+        
+    }
+}

src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.h

@@ -0,0 +1,12 @@
+#ifndef __SOLVE_H__
+#define __SOLVE_H__
+#include "Matrix.h"
+#include "Sparse.h"
+#include <vector>
+namespace Recon {
+    std::vector<std::vector<Aurora::Matrix>>
+    solveParameterIterator(Aurora::Sparse M, const Aurora::Matrix &b,
+                           const Aurora::Matrix &dims, bool oneIter = true, bool nonNeg = false);
+} // namespace Recon
+
+#endif // __SOLVE_H__