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

preprocessAScanBlockForReflection Test passed

Browse Source 
(left with Accuracy problem)
This commit is contained in:
kradchen
2023-06-13 14:24:05 +08:00
parent 9c98fcb899
commit 48f3d4d7f8
3 changed files with 506 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files

439
src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp Normal file
View File

@@ -0,0 +1,439 @@
#include "preprocessAScanBlockForReflection.h"
#include "Function.h"
#include "Function1D.h"
#include "Function2D.h"
#include "Function3D.h"
#include "Matrix.h"
#include "common/dataBlockCreation/removeDataFromArrays.h"
#include "config/config.h"
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstring>
#include <iostream>
#include <mkl_cblas.h>
#include <mkl_vml_functions.h>
#include <vector>
namespace Recon {
Aurora::Matrix _createDiffVector(const Aurora::Matrix & aMatrix, double beginValue,double endValue){
auto r = Aurora::zeros(1,aMatrix.getDataSize()+1);
r[0] = beginValue;
for (size_t i = 1; i < aMatrix.getDataSize(); i++)
{
r[i] = aMatrix[i] - aMatrix[i-1];
}
r[aMatrix.getDataSize()] = endValue;
return r;
}
/**
* @brief
*
* @param aMatrix
* @param op add 0, sub 1
* @return Aurora::Matrix
*/
Aurora::Matrix _blockOp(const Aurora::Matrix & aMatrix,int op){
if (aMatrix.isComplex()){
std::cerr<<"MatrixColDiff not support complex!"<<std::endl;
return Aurora::Matrix();
}
else{
if (aMatrix.isVector())
{
Aurora::Matrix result = Aurora::zeros(aMatrix.getDimSize(0)>1?aMatrix.getDataSize()-1:1,
aMatrix.getDimSize(0)>1?1:aMatrix.getDimSize(1)-1);
if (op == 1)vdSubI(result.getDataSize(),aMatrix.getData() , 1, aMatrix.getData()+1,1,result.getData(),1);
if (op == 0)vdAddI(result.getDataSize(),aMatrix.getData()+1 , 1, aMatrix.getData(),1,result.getData(),1);
return result;
}
else{
std::cerr<<"un expect operation!"<<std::endl;
return Aurora::Matrix();
// Aurora::Matrix result = Aurora::zeros(aMatrix.getDimSize(0)-1,aMatrix.getDimSize(1));
// #pragma omp parallel for
// for (size_t i = 0; i < aMatrix.getDimSize(1); i++)
// {
// if (op == 1){
// vdSubI(result.getDimSize(0),
// aMatrix.getData()+i*aMatrix.getDimSize(0)+1, 1,
// aMatrix.getData()+i*aMatrix.getDimSize(0)+1, 1,
// result.getData()+i*aMatrix.getDimSize(0),1);
// }
// if (op == 0){
// vdAddI(result.getDimSize(0),
// aMatrix.getData()+i*aMatrix.getDimSize(0)+1, 1,
// aMatrix.getData()+i*aMatrix.getDimSize(0)+1, 1,
// result.getData()+i*aMatrix.getDimSize(0),1);
// }
}
}
}
Aurora::Matrix generateGaussWindow(size_t winLength)
{
auto n = Aurora::linspace(0,5,winLength);
return exp( -0.2 * (transpose(n)^ 2));
}
//TODO:untested
Aurora::Matrix offsetSignalFourier(const PreComputes &preComputes, int nfft)
{
double dt = -(preComputes.offset/preComputes.timeInterval);
double binStart = floor((double)nfft/2);
auto temp = Aurora::zeros(1,nfft);
for (size_t i = 0; i < temp.getDataSize(); i++)
{
temp[i]=i;
}
temp= temp-binStart;
temp = Aurora::transpose(temp);
Aurora::ifftshift(temp);
auto fftBin = -1*dt*2*M_PI*temp/nfft;
auto negCfftbin = Aurora::complex(Aurora::zeros(fftBin.getDimSize(0),fftBin.getDimSize(1),fftBin.getDimSize(2)));
cblas_dcopy(fftBin.getDataSize(), fftBin.getData(), 1, negCfftbin.getData()+1, 2);
auto exponent = exp(negCfftbin);
if (!exponent.isComplex())exponent = complex(exponent);
return exponent;
}
std::vector<Aurora::Matrix> performSignalProcessing(
const Aurora::Matrix &blockedAScans,
const Aurora::Matrix &blockedSL,
const Aurora::Matrix &blockedRL,
const Aurora::Matrix &blockedSenderPosition,
const Aurora::Matrix &blockedReceiverPosition,
const Aurora::Matrix &blockedGain, const Aurora::Matrix &blockedChannels,
ExpInfo &info, const PreComputes &preComputes) {
std::vector<Aurora::Matrix> result;
auto t = size(blockedAScans);
size_t nSamples = (int)t[0];
size_t nAScans = (int)t[1];
auto valid = Aurora::ones(blockedAScans.getDimSize(1),1);
{
auto blockedGain_t = blockedGain.deepCopy();
Aurora::nantoval(blockedGain_t, 0);
valid = valid*blockedGain_t;
}
auto _blockedAScans = blockedAScans/blockedGain;
if (reflectParams::removeDCOffset == 1)
{
_blockedAScans = _blockedAScans-Aurora::mean(_blockedAScans);
}
int winLength = 100;
if (reflectParams::removeDCOffset == 1)
{
auto meanAS = mean(_blockedAScans.block(0,winLength-1,_blockedAScans.getDimSize(0) -1 -winLength));
_blockedAScans = _blockedAScans-meanAS.getScalar();
}
_blockedAScans.setBlockValue(0,0,winLength-1,0);
_blockedAScans.setBlockValue(0,_blockedAScans.getDimSize(0)-winLength-1,_blockedAScans.getDimSize(0)-1,0);
auto maxVal = max(abs(_blockedAScans));
auto stdVal = Aurora::std(_blockedAScans);
auto meanVal = Aurora::mean(abs(_blockedAScans));
auto ascanMapValue = meanVal*stdVal;
auto temp = (stdVal == 0) * (maxVal != 0);
Aurora::compareSet(valid,temp,1,0,Aurora::EQ);
if (reflectParams::findDefects==1)
{
valid = valid*(maxVal/meanVal >= reflectParams::epsilon);
}
// debugOutput.maxVal = maxVal;
// debugOutput.stdVal = stdVal;
// debugOutput.meanVal = meanVal;
// 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);
int count = reflectParams::suppressSameHeadLength* sizeof(double);
std::memset(begin_ptr,0,count);
}
}
}
auto fx = fft(_blockedAScans);
Aurora::Matrix fh;
if (reflectParams::useCorrelation && reflectParams::matchedFilterCeAScan)
{
if (!preComputes.measuredCEUsed)
{
fh = preComputes.matchedFilter.block(1, 0, nAScans-1);
}
else{
fh = preComputes.matchedFilter(Aurora::$,blockedChannels[0]-1).toMatrix();
}
double* value1 = Aurora::malloc(fx.getDataSize());
vdMulI(fx.getDataSize(), fx.getData(), 2, fh.getData(), 2, value1, 1);
double* value2 = Aurora::malloc(fx.getDataSize());
vdMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData() + 1, 2, value2, 1);
double* realData = Aurora::malloc(fx.getDataSize());
vdAdd(fx.getDataSize(), value1, value2, realData);
Aurora::Matrix real = Aurora::Matrix::New(realData, fx.getDimSize(0), fx.getDimSize(1));
vdMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData(), 2, value1, 1);
vdMulI(fx.getDataSize(), fx.getData(), 2, fh.getData() + 1, 2, value2, 1);
double* imagData = Aurora::malloc(fx.getDataSize());
vdSub(fx.getDataSize(), value1, value2, imagData);
Aurora::Matrix image = Aurora::Matrix::New(imagData, fx.getDimSize(0), fx.getDimSize(1));
double* complexData = Aurora::malloc(real.getDataSize(), true);
cblas_dcopy(real.getDataSize(), real.getData(), 1 , complexData ,2);
cblas_dcopy(image.getDataSize(), image.getData(), 1 , complexData + 1 ,2);
Aurora::Matrix complex = Aurora::Matrix::New(complexData, real.getDimSize(0), real.getDimSize(1), 1, Aurora::Complex);
_blockedAScans = Aurora::real(ifft(complex));
_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;
}
auto exponent = offsetSignalFourier(preComputes, nSamples);
// exponent.forceReshape(1, exponent.getDataSize(), 1);
_blockedAScans =_blockedAScans*exponent;
_blockedAScans = real(ifft(_blockedAScans));
if (reflectParams::removeTransmissionSignal == 1) {
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));
auto windowLength = reflectParams::windowLength;
auto start = round((distanceEPosRPos/reflectParams::expectedUSSpeedRange[1])*reflectParams::aScanReconstructionFrequency);
auto ende = round((distanceEPosRPos/reflectParams::expectedUSSpeedRange[0])*reflectParams::aScanReconstructionFrequency);
auto area = ceil(reflectParams::expectedPulseLength*info.Wavelength+1);
start = start- round(area*0.1);
ende = ende + round(area*0.9);
Aurora::compareSet(start, 1, 2, Aurora::NG);
Aurora::compareSet(start,nSamples,nSamples,Aurora::GT);
Aurora::compareSet(ende,nSamples,nSamples,Aurora::GT);
for (size_t idx = 0; idx < nAScans; idx++)
{
auto partData = _blockedAScans(Aurora::$,idx).toMatrix().block(0, start[idx]-1, ende[idx]-1);
auto l_partData = partData.getDataSize();
auto energyMove = _createDiffVector(partData,partData[0],0);
energyMove = _blockOp(energyMove,0);
energyMove = (energyMove^2)*0.5;
auto energyPot = partData^2;
auto energySum = energyMove + energyPot;
auto mean_energySum = Aurora::zeros(l_partData + windowLength,1);
for (size_t i = 0; i < windowLength; i++)
{
double* begin = mean_energySum.getData() + i;
int length = l_partData;
vdAddI(length, energySum.getData(), 1, begin, 1, begin, 1);
}
mean_energySum = mean_energySum.block(0,windowLength/2-1,mean_energySum.getDataSize()-1-round((windowLength-0.001)/2));
auto g = _createDiffVector(mean_energySum,0,0);
g = _blockOp(g,0);
auto f = _createDiffVector(g,g[0],0);
f = _blockOp(f,0);
Aurora::compareSet(f, 0, 0, Aurora::LT);
Aurora::compareSet(f,g,0,0,Aurora::NG);
Aurora::compareSet(f,0,0,Aurora::NG);
auto mean_energySum_f = f;
long scheitelRow=0,scheitelCol=0;
max(g,Aurora::Column,scheitelRow,scheitelCol);
long indexRow=0, indexCol=0;
//不知道对不对？？？
auto scheitel = scheitelRow+scheitelCol;
max(mean_energySum_f.block(0,0,scheitel),Aurora::Column,indexRow,indexCol);
size_t index = indexRow+indexCol;
double median_mean_energySum = Aurora::median(mean_energySum).getScalar();
bool flag = false;
if (mean_energySum[index] > median_mean_energySum)
{
for (size_t i = index; i > 0 ; i--)
{
if(mean_energySum[i]<median_mean_energySum){
index = i;
flag = true;
break;
}
}
if (!flag) index = 0;
}
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++)
{
if (hIndex_s[i]>0)hIndex.push_back(i);
}
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] ;
//在这里+1 调整会matlab的index便于后续计算的正确性
auto sig_end = index2 + start[idx] +1;
winLength = sig_end+winLength+200-sig_begin+1;
auto win = generateGaussWindow(winLength);
double* begin = _blockedAScans.getData() + (size_t)sig_begin-1;
int length = winLength;
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);
begin = _blockedAScans.getData() ;
length = sig_end;
vdMulI(length, begin, 1, win.getData(), 1, begin, 1);
}
}
_blockedAScans = fft(_blockedAScans);
if(reflectParams::useOptPulse==1){
auto n = nSamples;
auto nHalf = round((double)n/2);
_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 = abs(_blockedAScans.block(0,0,nSamples-1));
auto help = _blockedAScans.deepCopy();
help(0,Aurora::$) = 0;
help.setBlock(0, 1, help.getDimSize(0)-1, _blockOp(_blockedAScans, 1));
{
auto help_bbegin = help.block(0, 0, help.getDimSize(0)-2);
auto help_bend = help.block(0, 1, help.getDimSize(0)-1);
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;
}
help = _blockedAScans.deepCopy();
Aurora::compareSet(help,0,0,Aurora::LT);
if(reflectParams::limitNumPulsesTo>0)
{
if(size(help,1)>reflectParams::limitNumPulsesTo)
{
auto help2 = Aurora::sort(help);
auto temp = repmat(help2(help2.getDimSize(0)-reflectParams::limitNumPulsesTo,Aurora::$).toMatrix(),help.getDimSize(0),1);
Aurora::compareSet(help,temp,0,Aurora::LT);
}
}
if(reflectParams::normalizePeaks)
{
Aurora::compareSet(help,0,1,Aurora::GT);
}
help= real(ifft(fft(help)*(preComputes.sincPeak_ft.block(1,0,nAScans-1))));
Aurora::compareSet(_blockedAScans,0,help,Aurora::NL);
}
else{
_blockedAScans = real(ifft(_blockedAScans));
}
result.push_back(_blockedAScans);
result.push_back(valid);
ascanMapValue[0] = (float)ascanMapValue[0];
result.push_back(ascanMapValue);
return result;
}
preprocessAScanRResult preprocessAScanBlockForReflection(
Aurora::Matrix &AscanBlock, const Aurora::Matrix &blockedMP,
const Aurora::Matrix &blockedSL, const Aurora::Matrix &blockedSN,
const Aurora::Matrix &blockedRL, const Aurora::Matrix &blockedRN,
const Aurora::Matrix &blockedSenderPosition,
const Aurora::Matrix &blockedReceiverPosition,
const Aurora::Matrix &blockedGain, const Aurora::Matrix &blockedChannels,
ExpInfo &info, const PreComputes &preComputes)
{
preprocessAScanRResult result;
auto valid = Aurora::zeros(1, blockedMP.getDataSize()) ;
auto ascanMapValue = Aurora::zeros(1, blockedMP.getDataSize()) ;
// auto dOutMaxVal = Aurora::zeros(1, blockedMP.getDataSize()) ;
// auto dOutStdVal = Aurora::zeros(1, blockedMP.getDataSize()) ;
// auto dOutMeanVal = Aurora::zeros(1, blockedMP.getDataSize()) ;
// auto dOutSnr = Aurora::zeros(1, blockedMP.getDataSize()) ;
// auto dOutSnrPass = Aurora::zeros(1, blockedMP.getDataSize()) ;
switch (reflectParams::version) {
case 1:
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++)
{
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,
preComputes);
AscanBlock(Aurora::$,i) = signalPResult[0];
valid(Aurora::$,i) = signalPResult[1];
ascanMapValue(Aurora::$,i) = signalPResult[2];
// dOutMaxVal[idx] = dOut.maxVal;
// dOutStdVal[idx] = dOut.stdVal;
// dOutMeanVal[idx] = dOut.meanVal;
// dOutSnr[idx] = dOut.snr;
// dOutSnrPass[idx] = dOut.snrPass;
}
}
}
if (reflectParams::findDefects==1)
{
auto highNoiseScore = Aurora::zeros(1, ascanMapValue.getDataSize());
size_t count = 0;
for (size_t i = 0; i < valid.getDataSize(); i++)
{
if (valid[i])highNoiseScore[count++] = ascanMapValue[i];
}
highNoiseScore.forceReshape(1, count, 1);
double highNoiseScoreMean = mean(highNoiseScore).getScalar();
double highNoiseScoreStd = Aurora::std(highNoiseScore).getScalar();
double treshold99= (mean(highNoiseScore)+2.596*Aurora::std(highNoiseScore)).getScalar();
Aurora::compareSet(valid,ascanMapValue,treshold99,0,Aurora::GT);
// debugOutput.blockStatisticValue = ascanMapValue;
// debugOutput.blockStatisticUsedData = ascanMapValue > treshold99;
// debugOutput.maxVal = dOutMaxVal;
// debugOutput.stdVal = dOutStdVal;
// debugOutput.meanVal = dOutMeanVal;
// debugOutput.snr = dOutSnr;
// debugOutput.snrPass = dOutSnrPass;
}
result.AscanBlock = AscanBlock;
result.usedData = valid;
return result;
}
} // namespace Recon

26
src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.h Normal file
View File

@@ -0,0 +1,26 @@
#ifndef __PREPROCESSASCANBLOCKFORREFLECTION_H__
#define __PREPROCESSASCANBLOCKFORREFLECTION_H__
#include "Matrix.h"
#include "common/getMeasurementMetaData.h"
#include <cstddef>
namespace Recon {
struct preprocessAScanRResult {
Aurora::Matrix AscanBlock;
Aurora::Matrix usedData;
};
struct ExpInfo{
size_t Wavelength;
size_t expectedAScanLength;
};
preprocessAScanRResult preprocessAScanBlockForReflection(
Aurora::Matrix &AscanBlock, const Aurora::Matrix &mpBlock,
const Aurora::Matrix &slBlock, const Aurora::Matrix &snBlock,
const Aurora::Matrix &rlBlock, const Aurora::Matrix &rnBlock,
const Aurora::Matrix &senderPositionBlock,
const Aurora::Matrix &receiverPositionBlock,
const Aurora::Matrix &gainBlock, const Aurora::Matrix &channelBlock,
ExpInfo& info, const PreComputes &preComputes);
}
#endif // __PREPROCESSASCANBLOCKFORREFLECTION_H__

49
test/Reconstruction_Test.cpp
View File

@@ -8,6 +8,7 @@
#include "Matrix.h"
#include "Sparse.h"
#include "config/config.h"
#include "reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.h"
#include "transmissionReconstruction/reconstruction/buildMatrix/DGradient.h"
#include "transmissionReconstruction/reconstruction/buildMatrix/FMM.h"
#include "transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h"
@@ -46,16 +47,48 @@ protected:
};
TEST_F(Reconstruction_Test, determineOptimalPulse) {
Recon::reflectParams::imageResolution = 8.6381e-04;
Recon::reflectParams::optPulseFactor = 24;
auto result = Recon::determineOptimalPulse(1.0000e-07,3000);
EXPECT_EQ(3000, result.getDataSize());
ASSERT_DOUBLE_AE(result[2], 0.0025);
ASSERT_DOUBLE_AE(result[4], 0.0078);
Recon::reflectParams::imageResolution = 8.925316499483377e-04;
Recon::reflectParams::optPulseFactor = 48;
auto result = Recon::determineOptimalPulse(1.0000e-07,4096);
EXPECT_EQ(4096, result.getDataSize());
MatlabReader m2("/home/krad/TestData/sincPeakft.mat");
auto f1 = m2.read("sincPeak_ft");
for(size_t i=0; i<f1.getDataSize(); ++i)
{
EXPECT_DOUBLE_AE(f1[i], 0)<<"index:"<<i;
}
}
// for(size_t i=0; i<result.getDataSize(); ++i)
TEST_F(Reconstruction_Test, preprocessAScanBlockForReflection) {
Recon::initalizeConfig();
MatlabReader m2("/home/krad/TestData/preprocessRefC.mat");
auto blockedAScans = m2.read("blockedAScans");
auto blockedMP = m2.read("blockedMP");
auto blockedSL = m2.read("blockedSL");
auto blockedSN = m2.read("blockedSN");
auto blockedRL = m2.read("blockedRL");
auto blockedRN = m2.read("blockedRN");
auto senderPositionBlock = m2.read("blockedSenderPosition");
auto receiverPositionBlock = m2.read("blockedReceiverPosition");
auto channelBlock = m2.read("blockedChannels");
auto gainBlock = m2.read("blockedGain");
Recon::ExpInfo info;
info.expectedAScanLength= 4096;
info.Wavelength=3;
Recon::PreComputes preComputes;
preComputes.measuredCEUsed = true;
preComputes.offset = 5.2000e-07;
preComputes.timeInterval = 1.0000e-07;
preComputes.matchedFilter = m2.read("matchedFilter");
MatlabReader m("/home/krad/TestData/sincPeakft.mat");
preComputes.sincPeak_ft = m.read("sincPeak_ft");
auto result = Recon::preprocessAScanBlockForReflection(blockedAScans, blockedMP, blockedSL, blockedSN, blockedRL,
blockedRN,senderPositionBlock, receiverPositionBlock, gainBlock, channelBlock, info, preComputes);
size_t size = result.AscanBlock.getDataSize();
// for(size_t i=0; i<result.AscanBlock.getDataSize(); ++i)
// {
// ASSERT_DOUBLE_AE(f1[i], result.outSOS[i])<<"index:"<<i;
// ASSERT_DOUBLE_AE(f1[i], result[i])<<"index:"<<i;
// }
}