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

Improve preprocessAscanBlock

Browse Source
This commit is contained in:
kradchen
2023-06-19 09:33:38 +08:00
parent 420a1b5b1c
commit e9aaa0953d
1 changed files with 272 additions and 211 deletions
Download Patch File Download Diff File Expand all files Collapse all files

189
src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp
View File

@@ -12,6 +12,7 @@
#include <cstddef>
#include <cstring>
#include <iostream>
#include <memory>
#include <mkl_cblas.h>
#include <mkl_vml_functions.h>
#include <vector>
@@ -74,6 +75,12 @@ namespace Recon {
}
void fakeFree(void*){}
Aurora::Matrix getPartMatrixColRef(const Aurora::Matrix & aMatrix,int col){
std::shared_ptr<double> ptr{aMatrix.getData()+aMatrix.getDimSize(0)*col*aMatrix.getValueType(),fakeFree};
return Aurora::Matrix(ptr,{aMatrix.getDimSize(0),1,1},aMatrix.getValueType());
}
Aurora::Matrix generateGaussWindow(size_t winLength)
{
auto n = Aurora::linspace(0,5,winLength);
@@ -106,7 +113,7 @@ namespace Recon {
std::vector<Aurora::Matrix> performSignalProcessing(
const Aurora::Matrix &blockedAScans,
Aurora::Matrix &blockedAScans,
const Aurora::Matrix &blockedSL,
const Aurora::Matrix &blockedRL,
const Aurora::Matrix &blockedSenderPosition,
@@ -123,13 +130,18 @@ namespace Recon {
Aurora::nantoval(blockedGain_t, 0);
valid = valid*blockedGain_t;
}
auto _blockedAScans = blockedAScans/blockedGain;
int winLength = 100;
auto ascanMapValue = Aurora::zeros(1,blockedAScans.getDimSize(1),1);
#pragma omp parallel for num_threads(32)
for (size_t i = 0; i < blockedAScans.getDimSize(1); i++)
{
auto _blockedAScans = getPartMatrixColRef(blockedAScans,i);
_blockedAScans = _blockedAScans/blockedGain[i];
if (reflectParams::removeDCOffset == 1)
{
_blockedAScans = _blockedAScans-Aurora::mean(_blockedAScans);
_blockedAScans = _blockedAScans-Aurora::mean(_blockedAScans).getScalar();
}
int winLength = 100;
if (reflectParams::removeDCOffset == 1)
{
auto meanAS = mean(_blockedAScans.block(0,winLength-1,_blockedAScans.getDimSize(0) -1 -winLength));
@@ -140,12 +152,13 @@ namespace Recon {
auto maxVal = max(abs(_blockedAScans));
auto stdVal = Aurora::std(_blockedAScans);
auto meanVal = Aurora::mean(abs(_blockedAScans));
auto ascanMapValue = meanVal*stdVal;
ascanMapValue[i] = meanVal.getScalar()*stdVal.getScalar();
auto temp = (stdVal == 0) * (maxVal != 0);
Aurora::compareSet(valid,temp,1,0,Aurora::EQ);
valid[i] = temp[0]==1?0:valid[i];
// Aurora::compareSet(valid,temp,1,0,Aurora::EQ);
if (reflectParams::findDefects==1)
{
valid = valid*(maxVal/meanVal >= reflectParams::epsilon);
valid[i] = (valid[i]*(maxVal/meanVal >= reflectParams::epsilon)).getScalar();
}
// debugOutput.maxVal = maxVal;
// debugOutput.stdVal = stdVal;
@@ -153,28 +166,31 @@ namespace Recon {
// debugOutput.snr = maxVal ./ meanVal;
// debugOutput.snrPass = (maxVal./meanVal) < flags.dataSelection.epsilon;
if (reflectParams::suppressSameHead == 1)
{
for (size_t i = 0; i < blockedSL.getDataSize(); i++)
{
if (blockedSL[i] == blockedRL[i])
{
auto begin_ptr = _blockedAScans.getData()+i*_blockedAScans.getDimSize(0);
auto begin_ptr = _blockedAScans.getData();
int count = reflectParams::suppressSameHeadLength* sizeof(double);
std::memset(begin_ptr,0,count);
}
}
blockedAScans(Aurora::$,i)= _blockedAScans;
}
auto fx = fft(_blockedAScans);
Aurora::Matrix fh;
auto fx_all = fft(blockedAScans);
if (reflectParams::useCorrelation && reflectParams::matchedFilterCeAScan)
{
#pragma omp parallel for num_threads(32)
for (size_t i = 0; i < blockedAScans.getDimSize(1); i++)
{
Aurora::Matrix fx = getPartMatrixColRef(fx_all,i);
Aurora::Matrix fh;
if (!preComputes.measuredCEUsed)
{
fh = preComputes.matchedFilter.block(1, 0, nAScans-1);
}
else{
fh = preComputes.matchedFilter(Aurora::$,blockedChannels[0]-1).toMatrix();
fh = preComputes.matchedFilter(Aurora::$,blockedChannels[i]-1).toMatrix();
}
double* value1 = Aurora::malloc(fx.getDataSize());
vdMulI(fx.getDataSize(), fx.getData(), 2, fh.getData(), 2, value1, 1);
@@ -196,25 +212,36 @@ namespace Recon {
Aurora::Matrix complex = Aurora::Matrix::New(complexData, real.getDimSize(0), real.getDimSize(1), 1, Aurora::Complex);
Aurora::free(value1);
Aurora::free(value2);
_blockedAScans = Aurora::real(ifft(complex));
fx_all(Aurora::$,i) = complex;
}
blockedAScans = Aurora::real(ifft(fx_all));
_blockedAScans.setBlockValue(0, 0, winLength-1, 0);
_blockedAScans.setBlockValue(0, _blockedAScans.getDimSize(0)-winLength, _blockedAScans.getDimSize(0)-1, 0);
_blockedAScans = fft(_blockedAScans);
blockedAScans.setBlockValue(0, 0, winLength-1, 0);
blockedAScans.setBlockValue(0, blockedAScans.getDimSize(0)-winLength, blockedAScans.getDimSize(0)-1, 0);
blockedAScans = fft(blockedAScans);
}
else{
_blockedAScans = fx;
blockedAScans = fx_all;
}
auto exponent = offsetSignalFourier(preComputes, nSamples);
// exponent.forceReshape(1, exponent.getDataSize(), 1);
_blockedAScans =_blockedAScans*exponent;
_blockedAScans = real(ifft(_blockedAScans));
#pragma omp parallel for num_threads(32)
for (size_t i = 0; i < blockedAScans.getDimSize(1); i++)
{
Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i);
blockedAScans(Aurora::$,i) =_blockedAScans*exponent;
}
blockedAScans = real(ifft(blockedAScans));
if (reflectParams::removeTransmissionSignal == 1) {
#pragma omp parallel for num_threads(32)
for (size_t i = 0; i < blockedAScans.getDimSize(1); i++)
{
Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i);
auto distanceEPosRPos =
sqrt(
((blockedSenderPosition(0, Aurora::$).toMatrix() - blockedReceiverPosition(0, Aurora::$).toMatrix())^2) +
((blockedSenderPosition(1, Aurora::$).toMatrix() - blockedReceiverPosition(1, Aurora::$).toMatrix())^2) +
((blockedSenderPosition(2, Aurora::$).toMatrix() -blockedReceiverPosition(2, Aurora::$).toMatrix())^2));
((blockedSenderPosition(0, i).toMatrix() - blockedReceiverPosition(0, i).toMatrix())^2) +
((blockedSenderPosition(1, i).toMatrix() - blockedReceiverPosition(1, i).toMatrix())^2) +
((blockedSenderPosition(2, i).toMatrix() -blockedReceiverPosition(2, i).toMatrix())^2));
auto windowLength = reflectParams::windowLength;
auto start = round((distanceEPosRPos/reflectParams::expectedUSSpeedRange[1])*reflectParams::aScanReconstructionFrequency);
auto ende = round((distanceEPosRPos/reflectParams::expectedUSSpeedRange[0])*reflectParams::aScanReconstructionFrequency);
@@ -225,9 +252,9 @@ namespace Recon {
Aurora::compareSet(start,nSamples,nSamples,Aurora::GT);
Aurora::compareSet(ende,nSamples,nSamples,Aurora::GT);
for (size_t idx = 0; idx < nAScans; idx++)
// for (size_t idx = 0; idx < 1; idx++)
{
auto partData = _blockedAScans(Aurora::$,idx).toMatrix().block(0, start[idx]-1, ende[idx]-1);
auto partData = _blockedAScans.block(0, start[0]-1, ende[0]-1);
auto l_partData = partData.getDataSize();
auto energyMove = _createDiffVector(partData,partData[0],0);
@@ -238,9 +265,9 @@ namespace Recon {
auto energySum = energyMove + energyPot;
auto mean_energySum = Aurora::zeros(l_partData + windowLength,1);
for (size_t i = 0; i < windowLength; i++)
for (size_t k = 0; k < windowLength; k++)
{
double* begin = mean_energySum.getData() + i;
double* begin = mean_energySum.getData() + k;
int length = l_partData;
vdAddI(length, energySum.getData(), 1, begin, 1, begin, 1);
}
@@ -265,10 +292,10 @@ namespace Recon {
bool flag = false;
if (mean_energySum[index] > median_mean_energySum)
{
for (size_t i = index; i > 0 ; i--)
for (size_t k = index; k > 0 ; k--)
{
if(mean_energySum[i]<median_mean_energySum){
index = i;
if(mean_energySum[k]<median_mean_energySum){
index = k;
flag = true;
break;
}
@@ -277,46 +304,61 @@ namespace Recon {
}
auto hIndex_s =mean_energySum.block(0, scheitelCol, mean_energySum.getDataSize()-1)<(0.1*mean_energySum[scheitel]);
std::vector<int> hIndex;
for (size_t i = 0; i < hIndex_s.getDataSize(); i++)
for (size_t k = 0; k < hIndex_s.getDataSize(); k++)
{
if (hIndex_s[i]>0)hIndex.push_back(i);
if (hIndex_s[k]>0)hIndex.push_back(k);
}
if (hIndex.empty()){
hIndex.push_back(mean_energySum.getDataSize()-scheitel-1);
}
int index2 = hIndex[0]+scheitel-1+10;
//???存在没有index？？？matlab中有这一分支
auto sig_begin = index + start[idx] ;
auto sig_begin = index + start[0] ;
//在这里+1 调整会matlab的index便于后续计算的正确性
auto sig_end = index2 + start[idx] +1;
winLength = sig_end+winLength+200-sig_begin+1;
auto win = generateGaussWindow(winLength);
auto sig_end = index2 + start[0] +1;
int winLength2 = 100;
winLength2 = sig_end+winLength2+200-sig_begin+1;
auto win = generateGaussWindow(winLength2);
double* begin = _blockedAScans.getData() + (size_t)sig_begin-1;
int length = winLength;
int length = winLength2;
auto winReverse = Aurora::zeros(1,win.getDataSize());
std::reverse_copy(win.getData(), win.getData()+win.getDataSize(), winReverse.getData());
vdMulI(length, begin, 1, winReverse.getData(), 1, begin, 1);
winLength = sig_end;
win = generateGaussWindow(winLength);
winLength2 = sig_end;
win = generateGaussWindow(winLength2);
begin = _blockedAScans.getData() ;
length = sig_end;
vdMulI(length, begin, 1, win.getData(), 1, begin, 1);
// _blockedAScans = fft(_blockedAScans);
// blockedAScans(Aurora::$,i) =_blockedAScans;
}
}
_blockedAScans = fft(_blockedAScans);
}
// auto __blockedAScans = fft(blockedAScans(Aurora::$,0).toMatrix());
blockedAScans = fft(blockedAScans);
if(reflectParams::useOptPulse==1){
auto n = nSamples;
auto nHalf = round((double)n/2);
#pragma omp parallel for num_threads(32)
for (size_t i = 0; i < blockedAScans.getDimSize(1); i++)
{
Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i);
_blockedAScans(0,Aurora::$) = std::complex<double>{0,0};
_blockedAScans.setBlockComplexValue(0,nHalf,_blockedAScans.getDimSize(0)-1, std::complex<double>{0,0});
_blockedAScans.setBlock(0,1,nHalf-1,_blockedAScans.block(0,1,nHalf-1)*2);
}
blockedAScans = ifft(blockedAScans, n);
_blockedAScans = ifft(_blockedAScans, n);
blockedAScans = abs(blockedAScans);
_blockedAScans = abs(_blockedAScans.block(0,0,nSamples-1));
auto help = _blockedAScans.deepCopy();
auto help_all = blockedAScans.deepCopy();
#pragma omp parallel for num_threads(32)
for (size_t i = 0; i < blockedAScans.getDimSize(1); i++)
{
Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i);
Aurora::Matrix help = _blockedAScans.deepCopy();
help(0,Aurora::$) = 0;
help.setBlock(0, 1, help.getDimSize(0)-1, _blockOp(_blockedAScans, 1));
{
@@ -325,7 +367,7 @@ namespace Recon {
help_bend = (help_bend>0)*(help_bbegin<0);
auto tempBlock = Aurora::zeros(_blockedAScans.getDimSize(0),_blockedAScans.getDimSize(1),_blockedAScans.getDimSize(2));
tempBlock.setBlock(0, 0, tempBlock.getDimSize(0)-2, help_bend);
_blockedAScans = _blockedAScans*tempBlock;
blockedAScans(Aurora::$,i) = _blockedAScans*tempBlock;
}
help = _blockedAScans.deepCopy();
Aurora::compareSet(help,0,0,Aurora::LT);
@@ -347,20 +389,39 @@ namespace Recon {
{
Aurora::compareSet(help,0,1,Aurora::GT);
}
help_all(Aurora::$,i) = help;
}
help_all = fft(help_all);
#pragma omp parallel for num_threads(32)
for (size_t i = 0; i < help_all.getDimSize(1); i++)
{
Aurora::Matrix help = getPartMatrixColRef(help_all,i);
help_all(Aurora::$,i) = help*(preComputes.sincPeak_ft);
}
help= real(ifft(fft(help)*(preComputes.sincPeak_ft.block(1,0,nAScans-1))));
help_all = real(ifft(help_all));
#pragma omp parallel for num_threads(32)
for (size_t i = 0; i < help_all.getDimSize(1); i++)
{
Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i);
Aurora::Matrix help = getPartMatrixColRef(help_all,i);
Aurora::compareSet(_blockedAScans,0,help,Aurora::NL);
}
}
else{
_blockedAScans = real(ifft(_blockedAScans));
blockedAScans = real(ifft(blockedAScans));
}
result.push_back(_blockedAScans);
result.push_back(blockedAScans);
result.push_back(valid);
ascanMapValue[0] = (float)ascanMapValue[0];
for (size_t i = 0; i < ascanMapValue.getDataSize(); i++)
{
ascanMapValue[i] = (float)ascanMapValue[i];
}
result.push_back(ascanMapValue);
return result;
}
@@ -389,21 +450,21 @@ namespace Recon {
std::cerr<<"error USCT II only!!"<<std::endl;
break;
case 2:{
#pragma omp parallel for num_threads(32)
for (size_t i = 0; i < blockedMP.getDataSize(); i++)
// #pragma omp parallel for num_threads(32)
// for (size_t i = 0; i < blockedMP.getDataSize(); i++)
{
auto signalPResult = performSignalProcessing(
AscanBlock(Aurora::$, i).toMatrix(),
blockedSL(Aurora::$, i).toMatrix(),
blockedRL(Aurora::$, i).toMatrix(),
blockedSenderPosition(Aurora::$, i).toMatrix(),
blockedReceiverPosition(Aurora::$, i).toMatrix(),
blockedGain(Aurora::$, i).toMatrix(),
blockedChannels(Aurora::$, i).toMatrix(), info,
AscanBlock,
blockedSL,
blockedRL,
blockedSenderPosition,
blockedReceiverPosition,
blockedGain,
blockedChannels, info,
preComputes);
AscanBlock(Aurora::$,i) = signalPResult[0];
valid(Aurora::$,i) = signalPResult[1];
ascanMapValue(Aurora::$,i) = signalPResult[2];
AscanBlock = signalPResult[0];
valid = signalPResult[1];
ascanMapValue= signalPResult[2];
// dOutMaxVal[idx] = dOut.maxVal;
// dOutStdVal[idx] = dOut.stdVal;