feat: Add logic block log to performSignalProcessing in transmission reconstruction

This commit is contained in:
sunwen
2024-11-06 13:57:11 +08:00
parent b859728562
commit f4c16d121b
4 changed files with 43 additions and 23 deletions

View File

@@ -61,9 +61,9 @@ void producerThread( Parser* aParser, const Aurora::Matrix& aMotorPos,
int index=1;
for(int i=0; i<aMotorPos.getDataSize(); ++i)
{
for(int j=0; j<aSlList.getDataSize() / transParams::senderTASSize; ++j)
for(int j=0; j</*aSlList.getDataSize() / transParams::senderTASSize*/1; ++j)
{
for(int k=0; k<aSnList.getDataSize() / transParams::senderElementSize; ++k)
for(int k=0; k</*aSnList.getDataSize() / transParams::senderElementSize*/1; ++k)
{
// for(int i=0; i<1; ++i)
// {
@@ -74,8 +74,9 @@ void producerThread( Parser* aParser, const Aurora::Matrix& aMotorPos,
Matrix mp = aMotorPos(i).toMatrix();
Matrix sl = aSlList.block(0, transParams::senderTASSize*j, transParams::senderTASSize*j+transParams::senderTASSize - 1);
Matrix sn = aSnList.block(0, transParams::senderElementSize*k, transParams::senderElementSize*k+transParams::senderElementSize - 1);
RECON_INFO("getAscanBlockPreprocessed start");
auto blockData = getAscanBlockPreprocessed(aParser, mp, sl, sn, aRlList, aRnList, aGeom, aExpInfo, true, false);
RECON_INFO("getAscanBlockPreprocessed end");
float* channelListSizeData = Aurora::malloc(2);
channelListSizeData[0] = channelList.getDimSize(0);
channelListSizeData[1] = channelList.getDimSize(1);
@@ -88,9 +89,11 @@ void producerThread( Parser* aParser, const Aurora::Matrix& aMotorPos,
channelBlockData[i] = channelList[ind[i] - 1];
}
Matrix channelBlock = Matrix::New(channelBlockData, 1, channelBlockSize);
RECON_INFO("preprocessAScanBlockForReflection start");
auto preprocessData = preprocessAScanBlockForReflection(blockData.ascanBlockPreprocessed, blockData.mpBlock, blockData.slBlock,
blockData.snBlock, blockData.rlBlock, blockData.rnBlock, blockData.senderPositionBlock,
blockData.receiverPositionBlock, blockData.gainBlock, channelBlock, aExpInfo, aPreComputes);
RECON_INFO("preprocessAScanBlockForReflection end");
RECON_INFO("Reflect reconstructing...... {0}",index);
index++;
// std::unique_lock<std::mutex> lock(PRODUCER_MUTEX);

View File

@@ -13,6 +13,8 @@
#include "common/getMeasurementMetaData.h"
#include "config/config.h"
#include "log/log.h"
using namespace Aurora;
namespace Recon {
Matrix analyzeDetections(const Matrix &slBlockTotal, const Matrix &snBlockTotal, const Matrix &rlBlockTotal, const Matrix &rnBlockTotal, const Matrix &tofDataTotal){
@@ -154,6 +156,7 @@ namespace Recon {
startPosRef[i] = std::floor(std::max( tof2[i],(float)1.0));
endPosRef[i] = std::ceil(std::min(sizeAscan[1], tof2[i]+detectionWindowATT));
}
RECON_INFO("calculateAttenuation start");
return calculateAttenuation(envelopeOfAScan,startPos,endPos,envelopeOfReferenceAScan,startPosRef,endPosRef);
}
@@ -257,6 +260,7 @@ namespace Recon {
TimeWindowResult timeResult2;
timeResult2.AscanBlockProcessed = AscanRefBlock;
if (useTimeWindowing == 1) {
RECON_INFO("applyTimeWindowing start");
timeResult1 = applyTimeWindowing(
AscanBlock, sampleRate, distBlock, sosWaterBlock,
expectedSOSWater, offsetElectronicSamples, detectionWindowSOS,
@@ -267,6 +271,7 @@ namespace Recon {
minSpeedOfSound, maxSpeedOfSound, gaussWindow);
diffStartSearch = timeResult1.startSearch - timeResult2.startSearch;
RECON_INFO("applyTimeWindowing end");
}
auto _AscanBlock = timeResult1.AscanBlockProcessed;
auto _AscanRefBlock = timeResult2.AscanBlockProcessed;
@@ -277,12 +282,12 @@ namespace Recon {
auto mxl = std::min(maxlag, m - 1);
auto ceilLog2 = nextpow2(2 * m - 1);
auto m2 = pow(2, ceilLog2);
RECON_INFO("fft start");
auto x = fft(_AscanBlock, m2);
auto y = fft(_AscanRefBlock, m2);
auto c_1_1 = x * conj(y);
auto c_1_2 = ifft(c_1_1);
RECON_INFO("fft end");
auto c1 = real(c_1_2);
auto c = zeros(mxl + mxl + 1, c1.getDimSize(1));
c.setBlock(0, 0, mxl-1, c1.block(0, m2-mxl,m2-1));
@@ -309,11 +314,13 @@ namespace Recon {
DetectResult result;
result.sosValue = sosValue;
auto tofRel = tof - distBlock / sosWaterBlock;
RECON_INFO("detectAttVectorized start");
result.att = detectAttVectorized(
_AscanBlock, _AscanRefBlock, distRefBlock, sosWaterRefBlock,
tof, aScanReconstructionFrequency, offsetElectronic,
detectionWindowATT);
result.tof = tofRel;
RECON_INFO("detectAttVectorized end");
return result;
}

View File

@@ -57,12 +57,14 @@ namespace
std::mutex PROCESS_BUFFER_MUTEX;
std::condition_variable PROCESS_BUFFER_CONDITION;
int BUFFER_COUNT = 0;
int BUFFER_SIZE = 4;
int BUFFER_SIZE = 1;
Matrix prepareAScansForTransmissionDetection(const Matrix& aAscanBlock, const Matrix& aGainBlock)
{
RECON_INFO("prepareAScansForTransmissionDetection start");
Matrix result = aAscanBlock / repmat(aGainBlock, aAscanBlock.getDimSize(0), 1);
result = result - repmat(mean(result,FunctionDirection::Column), result.getDimSize(0), 1);
RECON_INFO("prepareAScansForTransmissionDetection end");
return result;
}
@@ -73,12 +75,17 @@ namespace
{
BlockOfTransmissionData result;
MetaInfos metaInfos;
RECON_INFO("getAscanBlockPreprocessed start");
auto blockData = getAscanBlockPreprocessed(aParser, aMp, aSl, aSn, aRlList, aRnList, aGeom, aExpInfo, true, true);
auto blockDataRef = getAscanBlockPreprocessed(aParserRef, aMpRef, aSl, aSn, aRlList, aRnList, aGeomRef, aExpInfoRef, true, true);
RECON_INFO("getAscanBlockPreprocessed end");
RECON_INFO("prepareAScansForTransmissionDetection start");
Matrix ascanBlock = prepareAScansForTransmissionDetection(blockData.ascanBlockPreprocessed, blockData.gainBlock);
Matrix ascanBlockRef = prepareAScansForTransmissionDetection(blockDataRef.ascanBlockPreprocessed, blockDataRef.gainBlock);
RECON_INFO("prepareAScansForTransmissionDetection end");
blockData.ascanBlockPreprocessed = Matrix();
blockDataRef.ascanBlockPreprocessed = Matrix();
RECON_INFO("Matrix * + start");
if(aExpInfo.Hardware == "USCT3dv3")
{
Matrix channelList = precalculateChannelList(aRlList, aRnList, aExpInfo, aPreComputes);
@@ -159,6 +166,7 @@ namespace
cblas_scopy(image.getDataSize(), image.getData(), 1 , complexData + 1 ,2);
complex = Matrix::New(complexData, real.getDimSize(0), real.getDimSize(1), 1, Aurora::Complex);
ascanBlockRef = Aurora::real(ifft(complex));
RECON_INFO("Matrix * + end");
}
else
{
@@ -168,16 +176,15 @@ namespace
if(transParams::applyCalib)
{
RECON_INFO("calculateSnr start");
metaInfos.snrValues = calculateSnr(ascanBlock, aSnrRmsNoise[0]);
metaInfos.snrValuesRef = calculateSnr(ascanBlockRef, aSnrRmsNoiseRef[0]);
RECON_INFO("calculateSnr end");
}
Matrix dists = distanceBetweenTwoPoints(blockData.senderPositionBlock, blockData.receiverPositionBlock);
Matrix distRefBlock = distanceBetweenTwoPoints(blockDataRef.senderPositionBlock, blockDataRef.receiverPositionBlock);
Matrix waterTempBlock = calculateWaterTemperature(aTasTemps.waterTempPreCalc_sl, aTasTemps.waterTempPreCalc_rl, blockData.slBlock, blockData.rlBlock, blockData.mpBlock);
Matrix waterTempRefBlock = calculateWaterTemperature(aTasTemps.waterTempRefPreCalc_sl, aTasTemps.waterTempRefPreCalc_rl, blockData.slBlock, blockData.rlBlock, blockDataRef.mpBlock);
if(transParams::saveDetection || transParams::outlierOnTasDetection || transParams::saveDebugInfomation)
{
metaInfos.mpBlock = blockData.mpBlock;
@@ -216,6 +223,7 @@ void getBlockOfTransmissionDataInThread(size_t aIndex, const Matrix& aMp, const
aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef);
std::unique_lock<std::mutex> lock(CREATE_BUFFER_MUTEX);
BLOCK_OF_TRANSIMISSIONDARA_BUFFER[std::to_string(aIndex)] = buffer;
RECON_INFO("add");
std::cout<<"add: "<<aIndex<<std::endl;
lock.unlock();
PROCESS_BUFFER_CONDITION.notify_one();
@@ -228,13 +236,13 @@ void createThreadForGetBlockOfTransmissionData(const Matrix& aMotorPos, const M
const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
{
size_t vectorSize = aMotorPos.getDataSize() * (aSlList.getDataSize() / transParams::senderTASSize) * (aSnList.getDataSize() / transParams::senderElementSize);
std::thread speedUpThread[vectorSize];
std::thread speedUpThread[/*vectorSize*/1];
for(int i=0; i<aMotorPos.getDataSize(); ++i)
{
for(int j=0; j<aSlList.getDataSize() / transParams::senderTASSize; ++j)
for(int j=0; j</*aSlList.getDataSize() / transParams::senderTASSize*/1; ++j)
{
for(int k=0; k<aSnList.getDataSize() / transParams::senderElementSize; ++k)
for(int k=0; k</*aSnList.getDataSize() / transParams::senderElementSize*/1; ++k)
{
// for(int i=0; i<1; ++i)
// {
@@ -252,8 +260,7 @@ void createThreadForGetBlockOfTransmissionData(const Matrix& aMotorPos, const M
std::unique_lock<std::mutex> lock(CREATE_BUFFER_MUTEX);
CREATE_BUFFER_CONDITION.wait(lock, []{return BUFFER_COUNT<BUFFER_SIZE;});
++BUFFER_COUNT;
lock.unlock();
std::cout<<"Remove: "<<index<<std::endl;
lock.unlock();
speedUpThread[index] = std::thread(getBlockOfTransmissionDataInThread,index,mp,mpRef,sl,sn,aRlList,aRnList,aTasTemps,aExpectedSOSWater,aGeom,aGeomRef,aSnrRmsNoise,aSnrRmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef);
}
}
@@ -278,12 +285,12 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
// % Load transmission detection data
// writeReconstructionLog('Loading transmission detection data. All available data from given motor positions are taken.', 1);
// [tofDataTotal, attDataTotal, senderList, receiverList, waterTempList, dataInfo] = loadTransmissionDetectionData(transParams.pathSaveDetection, transParams.pathData, motorPos, expInfo.rootMeasUniqueID);
TasTemps tasTemps;
tasTemps.waterTempPreCalc_rl = extractTasTemperature(aTemp.tasTemperature, aRlList, aMotorPos, aTemp.jumoTemp, transParams::minTemperature, transParams::maxTemperature);
tasTemps.waterTempPreCalc_sl = extractTasTemperature(aTemp.tasTemperature, aSlList, aMotorPos, aTemp.jumoTemp, transParams::minTemperature, transParams::maxTemperature);
tasTemps.waterTempRefPreCalc_rl = extractTasTemperature(aTempRef.tasTemperature, aRlList, aMotoPosRef, aTempRef.jumoTemp, transParams::minTemperature, transParams::maxTemperature);
tasTemps.waterTempRefPreCalc_sl = extractTasTemperature(aTempRef.tasTemperature, aSlList, aMotoPosRef, aTempRef.jumoTemp, transParams::minTemperature, transParams::maxTemperature);
aGeom.sensData = precalcSensitivity(aSlList, aSnList, aRlList, aRnList, aMotorPos, aGeom);
aGeomRef.sensData = aGeom.sensData;
@@ -293,7 +300,6 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
rmsNoise = estimateNoiseValueFromAScans(aSnList, aRnList, aGeom, aExpInfo, aParser);
rmsNoiseRef = estimateNoiseValueFromAScans(aSnList, aRnList, aGeom, aExpInfo, aParserRef);
}
size_t numScans = aMotorPos.getDataSize() * aSlList.getDataSize() * aSnList.getDataSize() * aRlList.getDataSize() * aRnList.getDataSize();
Matrix tofDataTotal = Matrix::fromRawData(new float[numScans], 1, numScans) + NAN;
Matrix attDataTotal = Matrix::fromRawData(new float[numScans], 1, numScans) + NAN;
@@ -316,7 +322,6 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
rlBlockTotal = zeros(1,numScans,1);
rnBlockTotal = zeros(1,numScans,1);
}
std::thread speedUpThread = std::thread(createThreadForGetBlockOfTransmissionData,aMotorPos,aMotoPosRef,aSlList,aSnList,aRlList,aRnList,tasTemps,aTemp.expectedSOSWater,aGeom,aGeomRef,rmsNoise,rmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef);
int numData = 0;
int numPossibleScans = 0;
@@ -331,9 +336,9 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
int finish_count=0;
for(int i=0; i<aMotorPos.getDataSize(); ++i)
{
for(int j=0; j<aSlList.getDataSize() / transParams::senderTASSize; ++j)
for(int j=0; j</*aSlList.getDataSize() / transParams::senderTASSize*/1; ++j)
{
for(int k=0; k<aSnList.getDataSize() / transParams::senderElementSize; ++k)
for(int k=0; k</*aSnList.getDataSize() / transParams::senderElementSize*/1; ++k)
{
// for(int i=0; i<1; ++i)
// {
@@ -349,14 +354,17 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
lock.unlock();
auto blockData = BLOCK_OF_TRANSIMISSIONDARA_BUFFER[std::to_string(index)];
RECON_INFO("transmissionDetection start");
DetectResult detect = transmissionDetection( blockData.ascanBlock, blockData.ascanBlockRef,
blockData.dists, blockData.distRefBlock,
blockData.waterTempBlock, blockData.waterTempRefBlock,
aTemp.expectedSOSWater[0]);
RECON_INFO("transmissionDetection end");
blockData.attData = detect.att;
blockData.tofData = detect.tof;
BlockOfTransmissionData transmissionBlock=blockData;
size_t numUsedData = transmissionBlock.senderBlock.getDimSize(1);
RECON_INFO("cblas_scopy start");
if(transParams::applyCalib)
{
cblas_scopy(numUsedData, transmissionBlock.metaInfos.snrValues.getData(), 1, snrValues.getData() + numData, 1);
@@ -378,14 +386,15 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
cblas_scopy(numUsedData, transmissionBlock.metaInfos.rlBlock.getData(), 1, rlBlockTotal.getData() + numData, 1);
cblas_scopy(numUsedData, transmissionBlock.metaInfos.rnBlock.getData(), 1, rnBlockTotal.getData() + numData, 1);
}
RECON_INFO("cblas_scopy end");
numData += numUsedData;
std::unique_lock<std::mutex> lockBufferCount(CREATE_BUFFER_MUTEX);
BLOCK_OF_TRANSIMISSIONDARA_BUFFER.erase(std::to_string(index));
--BUFFER_COUNT;
finish_count++;
lockBufferCount.unlock();
std::cout<<"\rTransmission reconstructing...... "<< finish_count<<"/64";
std::cout.flush();
RECON_INFO("Remove");
std::cout<<"Remove: "<<index<<std::endl;
CREATE_BUFFER_CONDITION.notify_one();
}
Recon::notifyProgress(6+10*((i+1)*(j+1)/(aMotorPos.getDataSize()*(aSlList.getDataSize()/ transParams::senderTASSize))));

View File

@@ -25,9 +25,10 @@ TransmissionReconstructionResult Recon::startTransmissionReconstruction(const Au
GeometryInfo& aGeomRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
{
RECON_INFO("Start getTransmissionData.");
RECON_INFO("getTransmissionData start");
auto transmissionData = getTransmissionData(aMotorPos, aMotoPosRef, aSlList, aSnList, aRlList, aRnList, aTemp, aTempRef,
aGeom, aGeomRef, aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef);
aGeom, aGeomRef, aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef);
RECON_INFO("getTransmissionData end");
Matrix dists = Recon::distanceBetweenTwoPoints(transmissionData.senderList, transmissionData.receiverList);
Matrix sosRef = Recon::temperatureToSoundSpeed(transmissionData.waterTempList, "marczak");
//Recon::notifyProgress(17);