From e57938ca5d8ef98f919c3817a446a4f0cc805087 Mon Sep 17 00:00:00 2001 From: kradchen Date: Wed, 14 Jun 2023 11:37:24 +0800 Subject: [PATCH] reconstructionSAFT 1 --- .../reconstructionSAFT/reconstructionSAFT.cpp | 93 +++++++------------ test/Reconstruction_Test.cpp | 46 +++++++++ 2 files changed, 82 insertions(+), 57 deletions(-) diff --git a/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp b/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp index 06c1c6f..90437a5 100644 --- a/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp +++ b/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp @@ -1,11 +1,13 @@ #include "reconstructionSAFT.h" #include "Function1D.h" +#include "Function2D.h" #include "Function3D.h" #include "Matrix.h" #include "SAFTStructs.h" #include "SAFT_ATT.h" #include "SAFT_TOFI.h" +#include "common/dataBlockCreation/removeDataFromArrays.h" #include "config/config.h" #include @@ -73,7 +75,7 @@ namespace Recon { Matrix_t receiverIdx{nullptr, 2, - {1, receiverIndex.getDataSize(), 1}, + { 1,receiverIndex.getDataSize(), 1}, receiverIndex.getDataSize()}; unsigned short *temp1 = new unsigned short[receiverIndex.getDataSize()]; std::copy(receiverIndex.getData(), @@ -83,7 +85,7 @@ namespace Recon { Matrix_t senderIdx{nullptr, 2, - {1, senderIndex.getDataSize(), 1}, + {senderIndex.getDataSize(),1, 1}, senderIndex.getDataSize()}; temp1 = new unsigned short[senderIndex.getDataSize()]; std::copy(senderIndex.getData(), @@ -118,8 +120,9 @@ namespace Recon { 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}; + int saftVariant_v[6]{(int)reflectParams::saftVariant[0],(int)reflectParams::saftVariant[1],(int)reflectParams::saftVariant[2], + (int)reflectParams::saftVariant[3],(int)reflectParams::saftVariant[4],(int)reflectParams::saftVariant[5]}; + Matrix_t saftVariant{saftVariant_v,2,{1,6,1},3}; params.push_back(saftVariant); Matrix_t SpeedMap3D{nullptr, @@ -197,9 +200,9 @@ namespace Recon { 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)}, + 2, + {1, (size_t)reconParams::gpuSelectionList.getDataSize(), + 1}, reconParams::gpuSelectionList.getDataSize()}; int *gdata = new int[reconParams::gpuSelectionList.getDataSize()]{0}; std::copy(reconParams::gpuSelectionList.getData(), reconParams::gpuSelectionList.getData() + reconParams::gpuSelectionList.getDataSize(), @@ -229,67 +232,43 @@ namespace Recon { recontructSAFT(const Aurora::Matrix &AScans, const Aurora::Matrix &senderPos, const Aurora::Matrix &receiverPos, const Aurora::Matrix &mpIndex, int SAFT_mode, TransRecos &transRecos, Aurora::Matrix &Env) { - auto TimeInterval = 1 / (reflectParams::aScanReconstructionFrequency); + double TimeInterval = 1.0 / (reflectParams::aScanReconstructionFrequency); std::vector motorPosAvailable; std::unique_copy(mpIndex.getData(), mpIndex.getData() + mpIndex.getDataSize(), std::back_inserter(motorPosAvailable)); for (auto mp : motorPosAvailable) { - std::vector 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 idsSender; - std::vector idsReceiver; - int addCount_s = 0; - int addCount_r = 0; + auto mpIdxs = mpIndex == mp; + size_t countMp = Aurora::sum(mpIdxs).getScalar(); - 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(); + auto senderPos_s = Aurora::transpose(removeDataFromArrays(senderPos, mpIdxs)); + auto receiverPos_s = Aurora::transpose(removeDataFromArrays(receiverPos, mpIdxs)); + Aurora::Matrix idsSender ; + auto senderPosGeom = Aurora::uniqueByRows(receiverPos_s, idsSender); + idsSender = Aurora::transpose(idsSender); + Aurora::Matrix idsReceiver ; + auto receiverPosGeom = Aurora::uniqueByRows(receiverPos_s, idsReceiver); + idsReceiver = Aurora::transpose(idsReceiver); + + auto numUsedData = countMp; 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[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, + 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); } } diff --git a/test/Reconstruction_Test.cpp b/test/Reconstruction_Test.cpp index 1f45dbc..c891cd9 100644 --- a/test/Reconstruction_Test.cpp +++ b/test/Reconstruction_Test.cpp @@ -3,6 +3,7 @@ #include #include "Function1D.h" +#include "Function3D.h" #include "MatlabReader.h" #include "MatlabWriter.h" #include "Matrix.h" @@ -16,6 +17,10 @@ #include "transmissionReconstruction/reconstruction/reconstruction.h" #include "transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h" +#include "reflectionReconstruction/preprocessData/preprocessTransmissionReconstructionForReflection.h" +#include "reflectionReconstruction/preprocessData/precalcImageParameters.h" +#include "reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.h" + #include "reflectionReconstruction/preprocessData/determineOptimalPulse.h" @@ -60,6 +65,47 @@ TEST_F(Reconstruction_Test, determineOptimalPulse) { } } +TEST_F(Reconstruction_Test, reconstructionSAFT) { + + Recon::initalizeConfig(); + MatlabReader m("/home/sun/testData/recontructSAFT.mat"); + auto Ascans = m.read("AScans"); + auto AttenuationMap3D = m.read("AttenuationMap3D"); + auto begin_TransMap = m.read("begin_TransMap"); + auto delta_TransMap = m.read("delta_TransMap"); + auto EnvOut = m.read("Env"); + auto mpIndex = m.read("mpIndex"); + auto receiverPos = m.read("receiverPos"); + auto SAFT_mode = m.read("SAFT_mode"); + auto senderPos = m.read("senderPos"); + auto SpeedMap3D = m.read("SpeedMap3D"); + Recon::TransRecos transRecos; + transRecos.speedMap3D = SpeedMap3D; + transRecos.attenuationMap3D = AttenuationMap3D; + transRecos.beginTransMap = begin_TransMap; + transRecos.deltaTransMap = delta_TransMap[0]; + MatlabReader m1("/home/sun/testData/resampleTransmissionVolume.mat"); + auto maxEmitter = m1.read("maxEmitter"); + auto maxReceiver = m1.read("maxReceiver"); + auto minEmitter = m1.read("minEmitter"); + auto minReceiver = m1.read("minReceiver"); + Recon::GeometryInfo geom; + geom.maxEmitter = maxEmitter; + geom.minEmitter = minEmitter; + geom.maxReceiver = maxReceiver; + geom.minReceiver = minReceiver; + geom.minSize = m1.read("minSize"); + geom.maxSize = m1.read("maxSize"); + precalcImageParameters(geom); + Aurora::Matrix Env = Aurora::zeros((int)Recon::reflectParams::imageXYZ[0],(int)Recon::reflectParams::imageXYZ[1],(int)Recon::reflectParams::imageXYZ[2]); + auto result = Recon::recontructSAFT(Ascans, senderPos, receiverPos, mpIndex, 2, transRecos, Env); + for(size_t i=0; i