Add preprocessTransmissionReconstructionForReflection.

This commit is contained in:
sunwen
2023-06-13 16:51:13 +08:00
parent dbb0353732
commit 15c90f9cef
4 changed files with 188 additions and 7 deletions

View File

@@ -112,7 +112,7 @@ namespace Recon {
const Aurora::Matrix &blockedSenderPosition, const Aurora::Matrix &blockedSenderPosition,
const Aurora::Matrix &blockedReceiverPosition, const Aurora::Matrix &blockedReceiverPosition,
const Aurora::Matrix &blockedGain, const Aurora::Matrix &blockedChannels, const Aurora::Matrix &blockedGain, const Aurora::Matrix &blockedChannels,
ExpInfo &info, const PreComputes &preComputes) { MeasurementInfo &info, const PreComputes &preComputes) {
std::vector<Aurora::Matrix> result; std::vector<Aurora::Matrix> result;
auto t = size(blockedAScans); auto t = size(blockedAScans);
size_t nSamples = (int)t[0]; size_t nSamples = (int)t[0];
@@ -371,7 +371,7 @@ namespace Recon {
const Aurora::Matrix &blockedSenderPosition, const Aurora::Matrix &blockedSenderPosition,
const Aurora::Matrix &blockedReceiverPosition, const Aurora::Matrix &blockedReceiverPosition,
const Aurora::Matrix &blockedGain, const Aurora::Matrix &blockedChannels, const Aurora::Matrix &blockedGain, const Aurora::Matrix &blockedChannels,
ExpInfo &info, const PreComputes &preComputes) MeasurementInfo &info, const PreComputes &preComputes)
{ {
preprocessAScanRResult result; preprocessAScanRResult result;
auto valid = Aurora::zeros(1, blockedMP.getDataSize()) ; auto valid = Aurora::zeros(1, blockedMP.getDataSize()) ;

View File

@@ -8,10 +8,6 @@ namespace Recon {
Aurora::Matrix AscanBlock; Aurora::Matrix AscanBlock;
Aurora::Matrix usedData; Aurora::Matrix usedData;
}; };
struct ExpInfo{
size_t Wavelength;
size_t expectedAScanLength;
};
preprocessAScanRResult preprocessAScanBlockForReflection( preprocessAScanRResult preprocessAScanBlockForReflection(
Aurora::Matrix &AscanBlock, const Aurora::Matrix &mpBlock, Aurora::Matrix &AscanBlock, const Aurora::Matrix &mpBlock,
@@ -20,7 +16,7 @@ namespace Recon {
const Aurora::Matrix &senderPositionBlock, const Aurora::Matrix &senderPositionBlock,
const Aurora::Matrix &receiverPositionBlock, const Aurora::Matrix &receiverPositionBlock,
const Aurora::Matrix &gainBlock, const Aurora::Matrix &channelBlock, const Aurora::Matrix &gainBlock, const Aurora::Matrix &channelBlock,
ExpInfo& info, const PreComputes &preComputes); MeasurementInfo& info, const PreComputes &preComputes);
} }
#endif // __PREPROCESSASCANBLOCKFORREFLECTION_H__ #endif // __PREPROCESSASCANBLOCKFORREFLECTION_H__

View File

@@ -0,0 +1,155 @@
#include "preprocessTransmissionReconstructionForReflection.h"
#include "Function3D.h"
#include "Matrix.h"
#include "reflectionReconstruction/preprocessData/resampleTransmissionVolume.h"
#include "src/config/config.h"
#include <algorithm>
#include <iostream>
#include <cstddef>
using namespace Aurora;
using namespace Recon;
PreprocessReflectionData Recon::preprocessTransmissionReconstructionForReflection(const Aurora::Matrix& aSOSMap, const Aurora::Matrix& aATTMap,
const Aurora::Matrix& aDdims, const GeometryInfo& aGeometryInfo,
const TempInfo& aTemp)
{
PreprocessReflectionData result;
int soundSpeedCorrection = reflectParams::soundSpeedCorrection;
int attenuationCorrection = reflectParams::attenuationCorrection;
Matrix minPos;
Matrix maxPos;
if(!aDdims.isNull())
{
minPos = Matrix::fromRawData(new double[3], 1, 3);
size_t maxPosSize = aDdims.getDataSize() - 3;
maxPos = Matrix::fromRawData(new double[maxPosSize], 1, maxPosSize);
std::copy(aDdims.getData(), aDdims.getData() + 3, minPos.getData());
std::copy(aDdims.getData() + 3, aDdims.getData() + 3 + maxPosSize, maxPos.getData());
}
if(soundSpeedCorrection == 1)
{
auto sosResampleVolume = resampleTransmissionVolume(aSOSMap, maxPos, minPos, aGeometryInfo);
result.transRecos.deltaTransMap = sosResampleVolume.deltaTransMap;
result.transRecos.beginTransMap = sosResampleVolume.beginTransMap;
// TemperatureModel4D暂无
// try % use temp model 4D
// speedmapx=transRecos.begin_TransMap(1):transRecos.delta_TransMap:(size(transRecos.SpeedMap3D,1)-1)*transRecos.delta_TransMap + transRecos.begin_TransMap(1);
// speedmapy=transRecos.begin_TransMap(2):transRecos.delta_TransMap:(size(transRecos.SpeedMap3D,2)-1)*transRecos.delta_TransMap + transRecos.begin_TransMap(2);
// speedmapz=transRecos.begin_TransMap(3):transRecos.delta_TransMap:(size(transRecos.SpeedMap3D,3)-1)*transRecos.delta_TransMap + transRecos.begin_TransMap(3);
// tidx=abs(transRecos.SpeedMap3D(:))<1; % voxels to fill
// [ix,iy,iz]=ind2sub(size(transRecos.SpeedMap3D),find(tidx));
// %fill with model
// if(~isempty(ix) && ~isempty(iy) && ~isempty(iz))
// transRecos.SpeedMap3D(tidx) = temperatureToSoundSpeed(polyvaln(temp.TemperatureModel4D.p,[speedmapx(ix)',speedmapy(iy)', speedmapz(iz)', zeros(size(ix))])); %surf(xg,yg,zg,reshape(tg,size(xg)))
// end
// catch
Matrix sppedMap3d = sosResampleVolume.transMap;
for(size_t i=0; i<sppedMap3d.getDataSize(); ++i)
{
if(std::abs(sppedMap3d[i]) < 1)
{
sppedMap3d[i] = aTemp.expectedSOSWater[0];
}
}
result.transRecos.speedMap3D = sppedMap3d;
//end
}
if( attenuationCorrection == 1)
{
Matrix attenuationMap3D = aATTMap.deepCopy();
for(size_t i=0; i<attenuationMap3D.getDataSize(); ++i)
{
if(attenuationMap3D[i] < 0)
{
attenuationMap3D[i] = 0;
}
}
attenuationMap3D = attenuationMap3D*20/std::log(10);
attenuationMap3D = attenuationMap3D*100;
auto attResampleVolume = resampleTransmissionVolume(attenuationMap3D, maxPos, minPos, aGeometryInfo);
result.transRecos.deltaTransMap = attResampleVolume.deltaTransMap;
result.transRecos.beginTransMap = attResampleVolume.beginTransMap;
attenuationMap3D = attResampleVolume.transMap;
if(result.transRecos.speedMap3D.isNull())
{
result.transRecos.speedMap3D = zeros(attenuationMap3D.getDimSize(0), attenuationMap3D.getDimSize(1), attenuationMap3D.getDimSize(2));
soundSpeedCorrection = 1;
}
result.transRecos.attenuationMap3D = attenuationMap3D;
}
else if(soundSpeedCorrection == 1 && !result.transRecos.speedMap3D.isNull())
{
Matrix speedMap3D = result.transRecos.speedMap3D;
result.transRecos.attenuationMap3D = zeros(speedMap3D.getDimSize(0), speedMap3D.getDimSize(1), speedMap3D.getDimSize(2));
attenuationCorrection = 1;
}
if(soundSpeedCorrection == 0 && attenuationCorrection == 0)
{
//与初始配置值不符暂不实现todo
// writeReconstructionLog('No sound speed correction / no attenuation corection -> using a water temperature based sound speed map and zero attenuation map.', 1);
// SAFT_mode = 2; % 0:(-SOS -ATT); 1:(+SOS -ATT); 2:(+SOS +ATT);
// delta_TransMap = flags.transmissionCorrection.resolutionTransmMap;
// minTransdPos = min(geometry.minEmitter,geometry.minReceiver);
// maxTransdPos =max(geometry.maxEmitter,geometry.maxReceiver);
// begin_TransMap = (min(flags.imageInfos.imageStartpoint',minTransdPos))' - 3*flags.transmissionCorrection.resolutionTransmMap; %elongate to cope for aperture movements
// end_SpeedMap = max(flags.imageInfos.imageEndpoint',maxTransdPos) + 3*flags.transmissionCorrection.resolutionTransmMap;
// end_SpeedMap(3) = end_SpeedMap(3) +0.025;
// x = (begin_TransMap(1):delta_TransMap:end_SpeedMap(1))+0.5*delta_TransMap;
// y = (begin_TransMap(2):delta_TransMap:end_SpeedMap(2))+0.5*delta_TransMap;
// z = (begin_TransMap(3):delta_TransMap:end_SpeedMap(3))+0.5*delta_TransMap;
// [XNEW, YNEW, ZNEW] = ndgrid(x,y,z);
// try
// transRecos.SpeedMap3D = reshape(temperatureToSoundSpeed(polyvaln(temp.TemperatureModel4D.p,[XNEW(:),YNEW(:), ZNEW(:), zeros(size(XNEW(:)))])),size(XNEW));
// catch
// try
// transRecos.SpeedMap3D = zeros([length(x), length(y), length(z)]) + temp.expectedSOSWater;
// catch
// writeReconstructionLog('No temperature data available. Can''t create an artificial sound speed map', 4);
// error([mfilename, ':noTemperatureData'],'No temperature data available. Can''t create an artificial sound speed map');
// end
// end
// transRecos.begin_TransMap = begin_TransMap;
// transRecos.delta_TransMap = delta_TransMap;
// transRecos.AttenuationMap3D = zeros(size(transRecos.SpeedMap3D));
// soundSpeedCorrection = 1; % setting sound speed correction to true for SAFT calls
}
else if (soundSpeedCorrection > 0 && attenuationCorrection == 0)
{
result.saftMode = 1;
result.transRecos.attenuationMap3D = zeros(result.transRecos.speedMap3D.getDimSize(0), result.transRecos.speedMap3D.getDimSize(1), result.transRecos.speedMap3D.getDimSize(2));
}
else if (soundSpeedCorrection > 0 && attenuationCorrection > 0)
{
if (result.transRecos.attenuationMap3D.getDimSize(0) != result.transRecos.speedMap3D.getDimSize(0) ||
result.transRecos.attenuationMap3D.getDimSize(1) != result.transRecos.speedMap3D.getDimSize(1) ||
result.transRecos.attenuationMap3D.getDimSize(2) != result.transRecos.speedMap3D.getDimSize(2))
{
std::cout<<"Size of SpeedMap3D ~= Size of AttenuationMap3D."<<std::endl;
}
result.saftMode = 2;
}
else if (soundSpeedCorrection == 0 && attenuationCorrection > 0)
{
std::cout<<"Only attenuation correction not allowed."<<std::endl;
}
result.soundSpeedCorrection = soundSpeedCorrection;
result.attenuationCorrection = attenuationCorrection;
return result;
}

View File

@@ -0,0 +1,30 @@
#ifndef PREPROCESS_TRANSMISSION_RECONSTRUCTION_FOR_REFLECTION
#define PREPROCESS_TRANSMISSION_RECONSTRUCTION_FOR_REFLECTION
#include "Matrix.h"
#include "src/common/getGeometryInfo.h"
#include "src/common/getMeasurementMetaData.h"
namespace Recon
{
struct TransRecos
{
Aurora::Matrix speedMap3D;
Aurora::Matrix attenuationMap3D;
Aurora::Matrix beginTransMap;
double deltaTransMap;
};
struct PreprocessReflectionData
{
TransRecos transRecos;
int saftMode;
double attenuationCorrection;
double soundSpeedCorrection;
};
PreprocessReflectionData preprocessTransmissionReconstructionForReflection(const Aurora::Matrix& aSOSMap, const Aurora::Matrix& aATTMap,
const Aurora::Matrix& aDdims, const GeometryInfo& aGeometryInfo,
const TempInfo& aTemp);
}
#endif