From c61a0ce09eaeca432233d496c77547e7459e6158 Mon Sep 17 00:00:00 2001
From: sunwen <sunwen@realcan.com>
Date: Thu, 1 Jun 2023 13:39:24 +0800
Subject: [PATCH] Update reconstructArt  and unittest. Fix buildMatrix bug.

---
 .../buildMatrix/buildMatrix.cpp               |  16 +--
 .../reconstruction/reconstruction.cpp         | 111 +++++++++++++++++-
 test/Reconstruction_Test.cpp                  |  45 +++++++
 3 files changed, 159 insertions(+), 13 deletions(-)

diff --git a/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp b/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp
index 9c32355..3a72f51 100644
--- a/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp
+++ b/src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp
@@ -1,5 +1,6 @@
 #include "buildMatrix.h"
 
+#include <cmath>
 #include <cstddef>
 #include <cstdio>
 #include <iostream>
@@ -34,12 +35,8 @@ namespace Recon
         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.path = round(path);
         result.weighting = getPixelLengthApproximation(startPt, endPt, res, result.pathLen);
         return result;
     }
@@ -304,19 +301,21 @@ namespace Recon
                 dims.getDimSize(1) == 2
                     ? convertToLinearIndices(dims, path.block(1, 0, 1))
                     : convertToLinearIndices(dims, path);
+            linearIndices = linearIndices-1;
+            linearIndices.forceReshape(1, linearIndices.getDataSize(), 1);
             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);
+            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);
+                i.setBlockValue(1, coeffIndex, cnt-1, rayCount-1);
+                j.setBlock(1,coeffIndex,cnt-1,linearIndices);
                 if (weighting.isScalar()){
                     s.setBlockValue(1, coeffIndex, cnt,weighting.getScalar());
                 }
@@ -343,6 +342,7 @@ namespace Recon
             }
             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;
     }
diff --git a/src/transmissionReconstruction/reconstruction/reconstruction.cpp b/src/transmissionReconstruction/reconstruction/reconstruction.cpp
index 3b29813..3fbdb05 100644
--- a/src/transmissionReconstruction/reconstruction/reconstruction.cpp
+++ b/src/transmissionReconstruction/reconstruction/reconstruction.cpp
@@ -7,10 +7,13 @@
 #include "config/config.h"
 
 #include "CudaEnvInit.h"
+#include "transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h"
+#include "src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.h"
 
 #include <algorithm>
 #include <cmath>
 #include <iostream>
+#include <vector>
 using namespace Aurora;
 namespace Recon {
     Aurora::Matrix calculateMinimalMaximalTransducerPositions(
@@ -115,19 +118,117 @@ namespace Recon {
                              Aurora::Matrix &receiverList, Aurora::Matrix &res,
                              double SOS_IN_WATER)
     {
-        auto nTotalRays = size(senderList, 2);
+        ArtResult result;
+        int nTotalRays = size(senderList, 2);
         int numIter = 1;
-        if (transParams::bentReconstruction)
+        bool bentRecon = transParams::bentReconstruction;
+        if (bentRecon)
         {
-            numIter =transParams::bentIter+1; 
+            numIter =transParams::bentIter + 1; 
         }
         for (size_t i = 0; i < transParams::gpuSelectionList.getDataSize(); i++)
         {
             std::string msg;
-            if (!resetGPUDevice((int)transParams::gpuSelectionList[i],msg))
+            // if (!resetGPUDevice((int)transParams::gpuSelectionList[i],msg))
+            // {
+            //     std::cerr<<msg<<std::endl;
+            // }
+        }
+        std::vector<int> potentialMapDataDims = {1,1,1};
+        for(size_t i=0; i<dims.getDataSize(); ++i)
+        {
+            potentialMapDataDims[i] = dims[i];
+        }
+        Matrix potentialMap = (zeros(potentialMapDataDims[0], potentialMapDataDims[1], potentialMapDataDims[2]) + 1) * SOS_IN_WATER;
+
+        Matrix b;
+        if(!data.isNull())
+        {
+            b = data;
+            b.forceReshape(b.getDataSize(), 1, 1);
+            for(size_t i=0; i<b.getDataSize(); ++i)
             {
-                std::cerr<<msg<<std::endl;
+                if(b[i] != b[i] || b[i] == INFINITY)
+                {
+                    b[i] = 0;
+                }
             }
         }
+        else if(bentRecon)
+        {
+            bentRecon = false;
+            numIter = 1;
+            std::cout<<"No SOS data. Attenuation reconstruction is carried out without bent ray calculations."<<std::endl;
+        }
+
+        Matrix bAtt;
+        if(!dataAtt.isNull())
+        {
+            bAtt = dataAtt;
+            bAtt.forceReshape(bAtt.getDataSize(), 1, 1);
+            for(size_t i=0; i<bAtt.getDataSize(); ++i)
+            {
+                if(bAtt[i] != bAtt[i] || bAtt[i] == INFINITY)
+                {
+                    bAtt[i] = 0;
+                }
+            }
+        }
+
+        std::vector<Matrix> allHitMaps;
+        if(transParams::saveDebugInfomation)
+        {
+            for(int i=0; i<numIter; ++i)
+            {
+                allHitMaps.push_back(Matrix());
+            }
+        }
+
+        transParams::nonNeg = false;
+        BuildMatrixResult buildMatrixR;
+        for(int iter=1; iter<=numIter; ++iter)
+        {
+            buildMatrixR = buildMatrix(senderList, receiverList, res, dims, bentRecon && (iter!=1), potentialMap);
+            if(!data.isNull() && bentRecon && iter != numIter)
+            {
+                //与默认配置bentRecon不符，暂不实现
+                // % reconstruction
+                // if(iter == 1)
+                //     solverOutput = solveParameterIterator(M, b, dims, parameters, 0); % straight ray with all mus and betas
+                // else
+                //     solverOutput = solveParameterIterator(M, b, dims, parameters, 1);
+                // end
+                // f1 = slownessToSOS(solverOutput{1}, SOS_IN_WATER);
+
+                // if iter == 1
+                //     outAll.straightRay = cellfun(@slownessToSOS, solverOutput, repmat({SOS_IN_WATER}, size(solverOutput)), 'UniformOutput', false);
+                // else
+                //     outAll.bentRay{iter-1} = f1;
+                // end
+                // % look if exit criterion reached
+                // if exitBent(f1, potentialMap, parameters, iter-1)
+                //     break;
+                // end
+
+                // potentialMap = f1;
+            }
+            else if(bentRecon && iter == numIter)
+            {
+                printf("Maximum of %d iterations for bent ray reconstruction reached. Exiting the loop.\n",iter);
+            }
+
+            if(transParams::saveDebugInfomation && !buildMatrixR.hitmap.isNull())
+            {
+                allHitMaps.push_back(buildMatrixR.hitmap);
+            }
+
+            if(!data.isNull())
+            {
+                result.outSOS = solveParameterIterator(buildMatrixR.M, b, dims, false, transParams::nonNeg)[0][0];
+            }
+        }
+
+        return result;
     }
+  
 } // namespace Recon
\ No newline at end of file
diff --git a/test/Reconstruction_Test.cpp b/test/Reconstruction_Test.cpp
index dbfdcf6..74c7fcf 100644
--- a/test/Reconstruction_Test.cpp
+++ b/test/Reconstruction_Test.cpp
@@ -39,6 +39,23 @@ protected:
     }
 };
 
+TEST_F(Reconstruction_Test, reconstructArt) {
+    MatlabReader m("/home/sun/testData/reconstructArt.mat");
+    auto data  = m.read("data");
+    auto dataAtt  = m.read("dataAtt");
+    auto dims = m.read("dims");
+    auto receiverList  = m.read("receiverList");
+    auto res = m.read("res");
+    auto senderList  = m.read("senderList");
+    auto SOS_IN_WATER  = m.read("SOS_IN_WATER");
+    auto f1 = m.read("f1");
+    auto result = Recon::reconstructArt(data, dataAtt, dims, senderList, receiverList, res, SOS_IN_WATER[0]);
+    for(size_t i=0; i<f1.getDataSize(); ++i)
+    {
+        EXPECT_DOUBLE_AE(f1[i], result.outSOS[i]);
+    }
+}
+
 TEST_F(Reconstruction_Test, checkAndScale) {
     MatlabReader m("/home/sun/testData/buildMatrix.mat");
     auto i  = m.read("i1");
@@ -51,6 +68,34 @@ TEST_F(Reconstruction_Test, checkAndScale) {
     Recon::checkAndScale(sparse,b,1196032);
 }
 
+TEST_F(Reconstruction_Test, buildMatrix) {
+    MatlabReader m("/home/sun/testData/buildMatrix.mat");
+    auto j1 = m.read("j1");
+    auto i1 = m.read("i1");
+    auto s1 = m.read("s1");
+    auto dims = m.read("dims");
+    auto receiverList  = m.read("receiverList");
+    auto potentialMap = m.read("potentialMap");
+    auto senderList  = m.read("senderList"); 
+    auto res  = m.read("res");
+    auto result = Recon::buildMatrix(senderList, receiverList, res, dims, false, potentialMap);
+    for(size_t i=0; i<i1.getDataSize(); ++i)
+    {
+        EXPECT_DOUBLE_AE(i1[i] - 1, result.M.getColVector()[i]);
+    }
+
+    for(size_t i=0; i<j1.getDataSize(); ++i)
+    {
+        EXPECT_DOUBLE_AE(j1[i]- 1, result.M.getRowVector()[i]);
+    }
+
+    for(size_t i=0; i<s1.getDataSize(); ++i)
+    {
+        EXPECT_DOUBLE_AE(s1[i], result.M.getValVector()[i]);
+    }
+
+}
+
 TEST_F(Reconstruction_Test, calculateMinimalMaximalTransducerPositions) {
     auto senderList  = Aurora::Matrix::fromRawData(new double[6]{1, 2, 3, 1, 2, 4}, 3, 2);
     auto receiverList  = Aurora::Matrix::fromRawData(new double[6]{1, 8, 3, 1, 2, 1}, 3, 2);