From 376b24597eb333a90c61ed2b4bfeaf9bc0913bb3 Mon Sep 17 00:00:00 2001 From: kradchen Date: Mon, 12 Jun 2023 16:53:19 +0800 Subject: [PATCH] Add start Refelection Reconstruction --- src/common/getMeasurementMetaData.h | 1 + .../reconstructionSAFT/reconstructionSAFT.cpp | 329 ++++++++++++++++++ .../reconstructionSAFT/reconstructionSAFT.h | 26 ++ .../startReflectionReconstruction.cpp | 60 ++++ .../startReflectionReconstruction.h | 14 + 5 files changed, 430 insertions(+) create mode 100644 src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp create mode 100644 src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.h create mode 100644 src/reflectionReconstruction/startReflectionReconstruction.cpp create mode 100644 src/reflectionReconstruction/startReflectionReconstruction.h diff --git a/src/common/getMeasurementMetaData.h b/src/common/getMeasurementMetaData.h index bf137f9..1961b30 100644 --- a/src/common/getMeasurementMetaData.h +++ b/src/common/getMeasurementMetaData.h @@ -64,6 +64,7 @@ namespace Recon bool measuredCEUsed; Aurora::Matrix measuredCE_TASIndices; Aurora::Matrix measuredCE_receiverIndices; + Aurora::Matrix sincPeak_ft; double offset; }; diff --git a/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp b/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp new file mode 100644 index 0000000..871d6c4 --- /dev/null +++ b/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp @@ -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 +#include +#include +#include +#include +#include + +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 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 motorPosAvailable; + std::unique_copy(mpIndex.getData(), mpIndex.getData() + mpIndex.getDataSize(), + motorPosAvailable.begin()); + 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; + + 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."< +#include +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<