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

Add start Refelection Reconstruction

Browse Source
This commit is contained in:
kradchen
2023-06-12 16:53:19 +08:00
parent 0c6b87fab3
commit 376b24597e
5 changed files with 430 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/common/getMeasurementMetaData.h
View File

@@ -64,6 +64,7 @@ namespace Recon
bool measuredCEUsed;
Aurora::Matrix measuredCE_TASIndices;
Aurora::Matrix measuredCE_receiverIndices;
Aurora::Matrix sincPeak_ft;
double offset;
};

329
src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp Normal file
View File

@@ -0,0 +1,329 @@
#include "reconstructionSAFT.h"
#include "Function1D.h"
#include "Function3D.h"
#include "Matrix.h"
#include "SAFTStructs.h"
#include "SAFT_ATT.h"
#include "SAFT_TOFI.h"
#include "config/config.h"
#include <algorithm>
#include <cstddef>
#include <cstdlib>
#include <iostream>
#include <sys/types.h>
#include <vector>
namespace Recon {
int compareColumn(const Aurora::Matrix &aMatrix, int colCount,
const Aurora::Matrix &aOther, int col2) {
for (size_t j = 0; j < colCount; j++) {
size_t i = 0;
for (; i < aMatrix.getDimSize(0); i++) {
if (aMatrix[j * aMatrix.getDimSize(0) + 1] !=
aOther[col2 * aOther.getDimSize(0) + i])
break;
}
if (i == aMatrix.getDimSize(0))
return j;
}
return -1;
}
Aurora::Matrix callSAFT(const Aurora::Matrix &AScans,
const Aurora::Matrix &receiverIndex,
const Aurora::Matrix &senderIndex,
const Aurora::Matrix &receiverPosGeom,
const Aurora::Matrix &senderPosGeom,
int SAFT_mode, double TimeInterval,
TransRecosParams transRecos,
const Aurora::Matrix &Env)
{
if (reflectParams::useAscanIndex == 0) {
if (reflectParams::soundSpeedCorrection == 0) {
std::cerr
<< "reconstruction without sound speed correction not allowed!!!!"
<< std::endl;
return Aurora::Matrix();
}
}
std::vector<Matrix_t> params;
Matrix_t AScan{nullptr,
(size_t)AScans.getDims(),
{(size_t)AScans.getDimSize(0), (size_t)AScans.getDimSize(1),
(size_t)AScans.getDimSize(2)},
AScans.getDataSize()};
float *ASdata = new float[AScans.getDataSize()]{0};
std::copy(AScans.getData(), AScans.getData() + AScans.getDataSize(),
ASdata);
AScan.Data = ASdata;
params.push_back(AScan);
Matrix_t imageStartpoint{nullptr, 2, {1, 3}, 3};
imageStartpoint.Data =
new float[3]{(float)reflectParams::imageStartpoint[0],
(float)reflectParams::imageStartpoint[1],
(float)reflectParams::imageStartpoint[2]};
params.push_back(imageStartpoint);
Matrix_t receiverIdx{nullptr,
2,
{1, receiverIndex.getDataSize(), 1},
receiverIndex.getDataSize()};
unsigned short *temp1 = new unsigned short[receiverIndex.getDataSize()];
std::copy(receiverIndex.getData(),
receiverIndex.getData() + receiverIndex.getDataSize(), temp1);
receiverIdx.Data = temp1;
params.push_back(receiverIdx);
Matrix_t senderIdx{nullptr,
2,
{1, senderIndex.getDataSize(), 1},
senderIndex.getDataSize()};
temp1 = new unsigned short[senderIndex.getDataSize()];
std::copy(senderIndex.getData(),
senderIndex.getData() + senderIndex.getDataSize(), temp1);
senderIdx.Data = temp1;
params.push_back(senderIdx);
Matrix_t receiverGeom{nullptr,
(size_t)receiverPosGeom.getDims(),
{(size_t)receiverPosGeom.getDimSize(0),
(size_t)receiverPosGeom.getDimSize(1),
(size_t)receiverPosGeom.getDimSize(2)},
receiverPosGeom.getDataSize()};
float *fdata = new float[receiverPosGeom.getDataSize()]{0};
std::copy(receiverPosGeom.getData(),
receiverPosGeom.getData() + receiverPosGeom.getDataSize(), fdata);
receiverGeom.Data = fdata;
params.push_back(receiverGeom);
Matrix_t senderGeom{nullptr,
(size_t)senderPosGeom.getDims(),
{(size_t)senderPosGeom.getDimSize(0),
(size_t)senderPosGeom.getDimSize(1),
(size_t)senderPosGeom.getDimSize(2)},
senderPosGeom.getDataSize()};
fdata = new float[senderPosGeom.getDataSize()]{0};
std::copy(senderPosGeom.getData(),
senderPosGeom.getData() + senderPosGeom.getDataSize(), fdata);
senderGeom.Data = fdata;
params.push_back(senderGeom);
Matrix_t _SAFT_mode{&SAFT_mode, 1, {1, 1, 1}, 1};
params.push_back(_SAFT_mode);
int saftVariant_v[3]{(int)reflectParams::saftVariant[0],(int)reflectParams::saftVariant[1],(int)reflectParams::saftVariant[2]};
Matrix_t saftVariant{saftVariant_v,2,{1,3,1},3};
params.push_back(saftVariant);
Matrix_t SpeedMap3D{nullptr,
(size_t)transRecos.SpeedMap3D.getDims(),
{(size_t)transRecos.SpeedMap3D.getDimSize(0),
(size_t)transRecos.SpeedMap3D.getDimSize(1),
(size_t)transRecos.SpeedMap3D.getDimSize(2)},
transRecos.SpeedMap3D.getDataSize()};
fdata = new float[transRecos.SpeedMap3D.getDataSize()]{0};
std::copy(transRecos.SpeedMap3D.getData(),
transRecos.SpeedMap3D.getData() + transRecos.SpeedMap3D.getDataSize(), fdata);
SpeedMap3D.Data = fdata;
params.push_back(SpeedMap3D);
Matrix_t BeginTransMap{nullptr,
(size_t)transRecos.BeginTransMap.getDims(),
{(size_t)transRecos.BeginTransMap.getDimSize(0),
(size_t)transRecos.BeginTransMap.getDimSize(1),
(size_t)transRecos.BeginTransMap.getDimSize(2)},
transRecos.BeginTransMap.getDataSize()};
fdata = new float[transRecos.BeginTransMap.getDataSize()]{0};
std::copy(transRecos.BeginTransMap.getData(),
transRecos.BeginTransMap.getData() + transRecos.BeginTransMap.getDataSize(), fdata);
BeginTransMap.Data = fdata;
params.push_back(BeginTransMap);
Matrix_t DeltaTransMap{nullptr,
(size_t)transRecos.DeltaTransMap.getDims(),
{(size_t)transRecos.DeltaTransMap.getDimSize(0),
(size_t)transRecos.DeltaTransMap.getDimSize(1),
(size_t)transRecos.DeltaTransMap.getDimSize(2)},
transRecos.DeltaTransMap.getDataSize()};
fdata = new float[transRecos.DeltaTransMap.getDataSize()]{0};
std::copy(transRecos.DeltaTransMap.getData(),
transRecos.DeltaTransMap.getData() + transRecos.DeltaTransMap.getDataSize(), fdata);
DeltaTransMap.Data = fdata;
params.push_back(DeltaTransMap);
Matrix_t AttenuationMap3D{nullptr,
(size_t)transRecos.AttenuationMap3D.getDims(),
{(size_t)transRecos.AttenuationMap3D.getDimSize(0),
(size_t)transRecos.AttenuationMap3D.getDimSize(1),
(size_t)transRecos.AttenuationMap3D.getDimSize(2)},
transRecos.AttenuationMap3D.getDataSize()};
fdata = new float[transRecos.AttenuationMap3D.getDataSize()]{0};
std::copy(transRecos.AttenuationMap3D.getData(),
transRecos.AttenuationMap3D.getData() + transRecos.AttenuationMap3D.getDataSize(), fdata);
AttenuationMap3D.Data = fdata;
params.push_back(AttenuationMap3D);
float imageResolution_v = (float)reflectParams::imageResolution;
Matrix_t imageResolution{&imageResolution_v, 1, {1, 1, 1}, 1};
params.push_back(imageResolution);
float TimeInterval_v = (float)TimeInterval;
Matrix_t _TimeInterval{&TimeInterval_v, 1, {1, 1, 1}, 1};
params.push_back(_TimeInterval);
int imageXYZ_v[3]{(int)reflectParams::imageXYZ[0],(int)reflectParams::imageXYZ[1],(int)reflectParams::imageXYZ[2]};
Matrix_t imageXYZ{imageXYZ_v,2,{1,3,1},3};
params.push_back(imageXYZ);
Matrix_t _Env{nullptr,
(size_t)Env.getDims(),
{(size_t)Env.getDimSize(0), (size_t)Env.getDimSize(1),
(size_t)Env.getDimSize(2)},
Env.getDataSize()};
float *Envdata = new float[Env.getDataSize()]{0};
std::copy(Env.getData(), Env.getData() + Env.getDataSize(),
Envdata);
_Env.Data = Envdata;
params.push_back(_Env);
int blockXYZ_v[3]{(int)reflectParams::blockDimXYZ[0],(int)reflectParams::blockDimXYZ[1],(int)reflectParams::blockDimXYZ[2]};
Matrix_t blockXYZ{blockXYZ_v,2,{1,3,1},3};
params.push_back(blockXYZ);
Matrix_t gpus{nullptr,
(size_t)reconParams::gpuSelectionList.getDims(),
{(size_t)reconParams::gpuSelectionList.getDimSize(0), (size_t)reconParams::gpuSelectionList.getDimSize(1),
(size_t)reconParams::gpuSelectionList.getDimSize(2)},
reconParams::gpuSelectionList.getDataSize()};
int *gdata = new int[reconParams::gpuSelectionList.getDataSize()]{0};
std::copy(reconParams::gpuSelectionList.getData(), reconParams::gpuSelectionList.getData() + reconParams::gpuSelectionList.getDataSize(),
gdata);
gpus.Data = gdata;
params.push_back(gpus);
if (reflectParams::useAscanIndex == 0){
Matrix_t medianWindowSize{&reflectParams::medianWindowSize,1,{1,1,1},1};
}
int other[2]={reflectParams::debugMode, reflectParams::attenuationCorrectionLimit} ;
Matrix_t otherP{other,2,{1,2,1},2};
params.push_back(otherP);
if (reflectParams::useAscanIndex == 0) {
auto result = SAFT_ATT(params);
return Aurora::Matrix::fromRawData((double *)result.Data, result.Dims[0],
result.Dims[1], result.Dims[2]);
} else {
auto result = SAFT_TOFI(params);
return Aurora::Matrix::fromRawData((double *)result.Data, result.Dims[0],
result.Dims[1], result.Dims[2]);
}
}
Aurora::Matrix
recontructSAFT(const Aurora::Matrix &AScans, const Aurora::Matrix &senderPos,
const Aurora::Matrix &receiverPos, const Aurora::Matrix &mpIndex,
int SAFT_mode, TransRecosParams &transRecos, Aurora::Matrix &Env) {
auto TimeInterval = 1 / (reflectParams::aScanReconstructionFrequency);
std::vector<int> motorPosAvailable;
std::unique_copy(mpIndex.getData(), mpIndex.getData() + mpIndex.getDataSize(),
motorPosAvailable.begin());
for (auto mp : motorPosAvailable) {
std::vector<int> mpIdxs;
for (size_t i = 0; i < mpIndex.getDataSize(); i++) {
if (mpIndex[i] == mp)
mpIdxs.push_back(i);
}
auto senderPosGeom = Aurora::zeros(senderPos.getDimSize(0), mpIdxs.size());
auto receiverPosGeom =
Aurora::zeros(receiverPos.getDimSize(0), mpIdxs.size());
std::vector<int> idsSender;
std::vector<int> idsReceiver;
int addCount_s = 0;
int addCount_r = 0;
for (auto mpIdx : mpIdxs) {
int flag = compareColumn(senderPosGeom, addCount_s, senderPos, mpIdx);
// old one, add ids
if (flag >= 0) {
idsSender.push_back(flag);
}
// found new unique one
else {
senderPosGeom(addCount_s, Aurora::$) = senderPos(mpIdx, Aurora::$);
idsSender.push_back(addCount_s++);
}
int flag2 =
compareColumn(receiverPosGeom, addCount_r, receiverPos, mpIdx);
// old one, add ids
if (flag2 >= 0) {
idsReceiver.push_back(flag2);
}
// found new unique one
else {
receiverPosGeom(addCount_r, Aurora::$) = receiverPos(mpIdx, Aurora::$);
idsReceiver.push_back(addCount_r++);
}
}
auto numUsedData = mpIdxs.size();
int count = numUsedData / reflectParams::blockSizeReco;
auto blockIdxs = Aurora::zeros(1, count + 2);
for (size_t i = 0; i < count + 1; i++) {
blockIdxs[i] = i * reflectParams::blockSizeReco;
}
blockIdxs[count + 1] = numUsedData;
for (size_t i = 0; i < count + 1; i++) {
size_t length = blockIdxs[i + 1] - blockIdxs[i];
size_t begin = blockIdxs[i];
size_t end = blockIdxs[i + 1];
auto _AScans = Aurora::zeros(AScans.getDimSize(0), length);
auto receiverIndex = Aurora::zeros(1, length);
auto senderIndex = Aurora::zeros(1, length);
for (size_t j = begin, k = 0; j < end; j++, k++) {
_AScans(Aurora::$, k) = AScans(Aurora::$, mpIdxs[j]);
receiverIndex[k] = idsReceiver[j];
senderIndex[k] = idsSender[j];
}
Env = callSAFT(_AScans, receiverIndex, senderIndex, receiverPosGeom,
senderPosGeom, SAFT_mode, TimeInterval, transRecos, Env);
}
}
return Env;
}
//TODO:untested
Aurora::Matrix polyvaln(PolyModel polymodel, const Aurora::Matrix &indepvar)
{
auto np = Aurora::size(indepvar);
int n = (int)np[0];
int p = (int)np[1];
Aurora::Matrix _indepvar = indepvar;
if (n == 1 && polymodel.ModelTerms.getDimSize(1) == 1){
_indepvar = Aurora::transpose(indepvar);
np = Aurora::size(indepvar);
n = (int)np[0];
p = (int)np[1];
}
else if(polymodel.ModelTerms.getDimSize(1) != p){
std::cerr<<"Size of indepvar array and this model are inconsistent."<<std::endl;
return Aurora::Matrix();
}
size_t nt = Aurora::size(polymodel.ModelTerms,1);
auto ypred = Aurora::zeros(n,1);
for (size_t i = 0; i < nt; i++)
{
auto t = Aurora::ones(n,1);
for (size_t j = 0; j < p; j++)
{
t = (t*(_indepvar(Aurora::$,j).toMatrix()^polymodel.Coefficients(i,j).toMatrix().getScalar()));
}
ypred = ypred+t*polymodel.Coefficients[i];
}
return ypred;
}
} // namespace Recon

26
src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.h Normal file
View File

@@ -0,0 +1,26 @@
#ifndef __RECONSTRUCTIONSAFT_H__
#define __RECONSTRUCTIONSAFT_H__
#include "Matrix.h"
namespace Recon {
struct TransRecosParams{
Aurora::Matrix SpeedMap3D;
Aurora::Matrix AttenuationMap3D;
Aurora::Matrix BeginTransMap;
Aurora::Matrix DeltaTransMap;
};
struct PolyModel{
Aurora::Matrix ModelTerms;
Aurora::Matrix Coefficients;
};
Aurora::Matrix
recontructSAFT(const Aurora::Matrix &AScans, const Aurora::Matrix &senderPos,
const Aurora::Matrix &receiverPos, const Aurora::Matrix &mpIndex,
int SAFT_mode, TransRecosParams &transRecos, Aurora::Matrix &Env);
Aurora::Matrix polyvaln(PolyModel polymodel, const Aurora::Matrix &indepvar);
}
#endif // __RECONSTRUCTIONSAFT_H__

60
src/reflectionReconstruction/startReflectionReconstruction.cpp Normal file
View File

@@ -0,0 +1,60 @@
#include "startReflectionReconstruction.h"
#include "Function3D.h"
#include "Matrix.h"
#include "common/precalculateChannelList.h"
#include "common/dataBlockCreation/getAScanBlockPreprocessed.h"
#include "reflectionReconstruction/preprocessData/determineOptimalPulse.h"
#include "config/config.h"
#include "CudaEnvInit.h"
#include <cstdio>
#include <iostream>
namespace Recon
{
Aurora::Matrix startReflectionReconstruction(Parser* aParser, int SAFT_mode, const Aurora::Matrix &motorPos,
const Aurora::Matrix &slList, const Aurora::Matrix &snList,
const Aurora::Matrix &rlList, const Aurora::Matrix &rnList,
const Aurora::Matrix &geom, MeasurementInfo expInfo, PreComputes preComputes)
{
printf("Reflection reconstruction is carried out.");
printf("Preperations for reconstructions.");
printf(" - reset GPUs");
for (size_t i = 0; i < reconParams::gpuSelectionList.getDataSize(); i++)
{
std::string msg;
if (!resetGPUDevice((int)transParams::gpuSelectionList[i],msg))
{
std::cerr<<msg<<std::endl;
}
}
printf(" - optimal pulse");
if(reflectParams::useOptPulse==1 && reflectParams::runReflectionReco){
preComputes.sincPeak_ft = determineOptimalPulse(preComputes.timeInterval ,expInfo.expectedAScanLength);
}
printf(" - channel list");
auto channelList = precalculateChannelList(rlList, rnList, expInfo, preComputes);
size_t numScans = motorPos.getDataSize() * slList.getDataSize() *
snList.getDataSize() * rlList.getDataSize() *
rnList.getDataSize();
printf(" - blocking");
auto Env = Aurora::zeros((int)reflectParams::imageXYZ[0],(int)reflectParams::imageXYZ[1],(int)reflectParams::imageXYZ[2]);
int numTakenScans = 0,numProcessedScans = 0,numPossibleScans = 0;
for(int i=0; i<motorPos.getDataSize(); ++i)
{
#pragma omp parallel for num_threads(24)
for(int j=0; j<slList.getDataSize() / transParams::senderTASSize; ++j)
{
for(int k=0; k<snList.getDataSize() / transParams::senderElementSize; ++k)
{
Aurora::Matrix mp = motorPos(i).toMatrix();
}
}
}
return Aurora::Matrix();
}
}

14
src/reflectionReconstruction/startReflectionReconstruction.h Normal file
View File

@@ -0,0 +1,14 @@
#ifndef __STARTREFLECTIONRECONSTRUCTION_H__
#define __STARTREFLECTIONRECONSTRUCTION_H__
#include "Matrix.h"
#include "transmissionReconstruction/detection/getTransmissionData.h"
#include "reconstructionSAFT/reconstructionSAFT.h"
namespace Recon {
Aurora::Matrix startReflectionReconstruction(
Parser* aParser, int SAFT_mode, const Aurora::Matrix &motorPos,
const Aurora::Matrix &slList, const Aurora::Matrix &snList,
const Aurora::Matrix &rlList, const Aurora::Matrix &rnList,
const Aurora::Matrix &geom,TransRecosParams &transRecos,
MeasurementInfo expInfo, PreComputes preComputes);
}
#endif // __STARTREFLECTIONRECONSTRUCTION_H__