234 lines
8.5 KiB
C++
234 lines
8.5 KiB
C++
#include "reconstruction.h"
|
||
#include "Function.h"
|
||
#include "Function1D.h"
|
||
#include "Function2D.h"
|
||
#include "Function3D.h"
|
||
#include "Matrix.h"
|
||
#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(
|
||
const Aurora::Matrix &aMSenderList, const Aurora::Matrix &aMReceiverList) {
|
||
|
||
auto minEmitter = min(aMSenderList, Row);
|
||
auto maxEmitter = max(aMSenderList, Row);
|
||
|
||
auto minReceiver = min(aMReceiverList, Row);
|
||
auto maxReceiver = max(aMReceiverList, Row);
|
||
|
||
auto data1 = Aurora::malloc(minEmitter.getDataSize() * 2);
|
||
auto data2 = Aurora::malloc(minEmitter.getDataSize() * 2);
|
||
|
||
auto minM = Matrix::New(data1,minEmitter.getDimSize(0),2);
|
||
auto maxM = Matrix::New(data2,minEmitter.getDimSize(0),2);
|
||
minM($,0) = minEmitter;
|
||
minM($,1) = minReceiver;
|
||
maxM($,0) = maxEmitter;
|
||
maxM($,1) = maxReceiver;
|
||
minM = min(minM,Row);
|
||
maxM = max(maxM,Row);
|
||
auto data3 = Aurora::malloc(minM.getDataSize() * 2);
|
||
auto ddims = Matrix::New(data3,minM.getDimSize(0),2);
|
||
ddims($,0) = minM;
|
||
ddims($,1) = maxM;
|
||
ddims.forceReshape(1, ddims.getDataSize(), 1);
|
||
return ddims;
|
||
}
|
||
|
||
Aurora::Matrix calculateResolution(const Aurora::Matrix &aVDdims, const Aurora::Matrix &aVDims)
|
||
{
|
||
auto numDim = aVDims.getDataSize();
|
||
auto res = (aVDdims.block(1,numDim,aVDdims.getDataSize()-1) - aVDdims.block(1,0,numDim-1))/ (aVDims - 1);
|
||
|
||
if(numDim == 3 && aVDims[2] == 1){
|
||
res[2] = res[0];
|
||
}
|
||
return res;
|
||
}
|
||
|
||
Aurora::Matrix getDimensions(double aNumPixelXY, const Aurora::Matrix& aVDdims)
|
||
{
|
||
int numDim = aVDdims.getDataSize()/2;
|
||
|
||
if (!(numDim == 2 || numDim == 3)){
|
||
std::cerr<<"Inputs does not match requirements."<<std::endl;
|
||
return Aurora::Matrix();
|
||
}
|
||
|
||
auto dims = ones(1,numDim);
|
||
|
||
dims[0] = aNumPixelXY;
|
||
dims[1] = aNumPixelXY;
|
||
|
||
if (numDim == 3)
|
||
{
|
||
dims[2] = std::ceil((aVDdims[5]-aVDdims[2])/((aVDdims[3]-aVDdims[0])/dims[0]));
|
||
//BUG:发现源代码可能存在无限值，加了std::isfinite验证
|
||
if(dims[2] == 0 || !std::isfinite(dims[2]))dims[2] = 1;
|
||
}
|
||
return dims;
|
||
}
|
||
|
||
void slownessToSOS(Aurora::Matrix & aVF1, double aSOS_IN_WATER)
|
||
{
|
||
aVF1 = 1 / ((aVF1 + 1) / aSOS_IN_WATER);
|
||
}
|
||
|
||
DiscretizePositionValues discretizePositions(Aurora::Matrix &aVSenderCoordList, Aurora::Matrix &aVReceiverCoordList, double aNumPixelXY)
|
||
{
|
||
DiscretizePositionValues result;
|
||
result.ddims = calculateMinimalMaximalTransducerPositions(aVSenderCoordList, aVReceiverCoordList);
|
||
result.dims = getDimensions(aNumPixelXY, result.ddims);
|
||
result.res = calculateResolution(result.ddims, result.dims);
|
||
int numDim = aVSenderCoordList.getDimSize(0);
|
||
auto tempDdims = result.ddims.block(0,0,numDim-1);
|
||
auto tempTransRes = repmat(transpose(result.res),1,aVSenderCoordList.getDimSize(1));
|
||
auto tempDivRes = repmat(transpose(tempDdims / result.res),1,aVSenderCoordList.getDimSize(1));
|
||
aVSenderCoordList = (aVSenderCoordList / tempTransRes) - tempDivRes + 1;
|
||
aVReceiverCoordList = (aVReceiverCoordList / tempTransRes) - tempDivRes + 1;
|
||
{
|
||
float* temp = new float[aVReceiverCoordList.getDataSize()]{0};
|
||
std::copy(aVReceiverCoordList.getData(),aVReceiverCoordList.getData()+aVReceiverCoordList.getDataSize(),temp);
|
||
std::copy(temp,temp+aVReceiverCoordList.getDataSize(),aVReceiverCoordList.getData());
|
||
delete [] temp;
|
||
}
|
||
{
|
||
float* temp = new float[aVSenderCoordList.getDataSize()]{0};
|
||
std::copy(aVSenderCoordList.getData(),aVSenderCoordList.getData()+aVSenderCoordList.getDataSize(),temp);
|
||
std::copy(temp,temp+aVSenderCoordList.getDataSize(),aVSenderCoordList.getData());
|
||
delete [] temp;
|
||
}
|
||
result.receiverCoordList = aVReceiverCoordList;
|
||
result.senderCoordList = aVSenderCoordList;
|
||
return result;
|
||
}
|
||
|
||
ArtResult reconstructArt(Aurora::Matrix &data, Aurora::Matrix &dataAtt,
|
||
Aurora::Matrix &dims,
|
||
Aurora::Matrix &senderList,
|
||
Aurora::Matrix &receiverList, Aurora::Matrix &res,
|
||
double SOS_IN_WATER)
|
||
{
|
||
ArtResult result;
|
||
int nTotalRays = size(senderList, 2);
|
||
int numIter = 1;
|
||
bool bentRecon = transParams::bentReconstruction;
|
||
if (bentRecon)
|
||
{
|
||
numIter =transParams::bentIter + 1;
|
||
}
|
||
for (size_t i = 0; i < transParams::gpuSelectionList.getDataSize(); i++)
|
||
{
|
||
std::string 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)
|
||
{
|
||
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
|