sunwen/UR
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Update reconstructArt and unittest.

Browse Source 
Fix buildMatrix bug.
This commit is contained in:
sunwen
2023-06-01 13:39:24 +08:00
parent 367176f86c
commit c61a0ce09e
3 changed files with 159 additions and 13 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
src/transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cpp
View File

@@ -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;
}

111
src/transmissionReconstruction/reconstruction/reconstruction.cpp
View File

@@ -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

45
test/Reconstruction_Test.cpp
View File

@@ -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);