First update.

This commit is contained in:
sunwen
2023-10-09 09:29:21 +08:00
parent a434c105c0
commit 903ac2c087
56 changed files with 1044 additions and 1021 deletions

View File

@@ -10,9 +10,9 @@
#include <omp.h>
namespace Recon {
Aurora::Matrix determineOptimalPulse(double timeInterval, size_t expectedAScanLength)
Aurora::Matrix determineOptimalPulse(float timeInterval, size_t expectedAScanLength)
{
double imgResolution = reflectParams::imageResolution;
float imgResolution = reflectParams::imageResolution;
int optPulseFactor = reflectParams::optPulseFactor;
auto minOptPulse = ceil(8*imgResolution*(1/timeInterval)/1500);
@@ -27,11 +27,11 @@ namespace Recon {
}
auto desSamplingFreq = (1/timeInterval)*0.5;
double sincFact = optPulseFactor;
double sincLength = round(11*sincFact/16);
float sincFact = optPulseFactor;
float sincLength = round(11*sincFact/16);
sincLength = 3*sincLength;
double begin =-sincLength/2+(timeInterval)/2;
double end = sincLength/2;
float begin =-sincLength/2+(timeInterval)/2;
float end = sincLength/2;
int length = end - begin + 1;
auto sincT = Aurora::zeros(1,length);
for (int j = 0; begin <= end; begin++, j++)
@@ -47,15 +47,15 @@ namespace Recon {
Aurora::padding(sincPeak, expectedAScanLength-1, 0);
// sincPeak = Aurora::transpose(sincPeak);
size_t offset = floor((double)sincPeak_len/2);
size_t offset = floor((float)sincPeak_len/2);
//cicshift
// #pragma omp parallel for
for (size_t i = 0; i < sincPeak.getDimSize(1); i++)
{
double *temp = new double[offset]{0};
double * beginPtr = sincPeak.getData()+i*sincPeak.getDimSize(0);
double * endPtr = beginPtr + offset;
double * col_end = beginPtr + (i+1)*sincPeak.getDimSize(0);
float *temp = new float[offset]{0};
float * beginPtr = sincPeak.getData()+i*sincPeak.getDimSize(0);
float * endPtr = beginPtr + offset;
float * col_end = beginPtr + (i+1)*sincPeak.getDimSize(0);
std::copy(beginPtr,endPtr,temp);
std::copy(endPtr,col_end,beginPtr);
std::copy(temp,temp+offset,col_end-offset);

View File

@@ -3,7 +3,7 @@
#include "Matrix.h"
namespace Recon {
Aurora::Matrix determineOptimalPulse(double timeInterval,
Aurora::Matrix determineOptimalPulse(float timeInterval,
size_t expectedAScanLength);
}

View File

@@ -10,11 +10,11 @@
using namespace Recon;
using namespace Aurora;
double Recon::estimateOffset(const MeasurementInfo aExpInfo, const CEInfo& aCe, const Aurora::Matrix& aMatchedFilter)
float Recon::estimateOffset(const MeasurementInfo aExpInfo, const CEInfo& aCe, const Aurora::Matrix& aMatchedFilter)
{
double offset = 0;
float offset = 0;
Matrix ce = aCe.ceRef;
double ceOffSet = aCe.ceRefOffSet;
float ceOffSet = aCe.ceRefOffSet;
Matrix ceMeasured;
if(aCe.measuredCEAvailable)
{
@@ -35,7 +35,7 @@ double Recon::estimateOffset(const MeasurementInfo aExpInfo, const CEInfo& aCe,
{
size_t column = ceMeasured.getDimSize(1);
size_t ceRows = ce.getDimSize(0);
double* offsetCEMeasuredData = Aurora::malloc(column);
float* offsetCEMeasuredData = Aurora::malloc(column);
for(int i=0; i<column; ++i)
{
Matrix corrCE = xcorr(ceMeasured($,i).toMatrix(), ce);
@@ -54,9 +54,9 @@ double Recon::estimateOffset(const MeasurementInfo aExpInfo, const CEInfo& aCe,
}
else
{
double dacDelay = aExpInfo.dacDelay;
double filterDisabled = aExpInfo.filterByPass;
double dacDelayInS = dacDelay / 10000000;
float dacDelay = aExpInfo.dacDelay;
float filterDisabled = aExpInfo.filterByPass;
float dacDelayInS = dacDelay / 10000000;
if(filterDisabled == 1)
{
offset += aCe.offsetFilterDisabled + dacDelayInS;
@@ -76,7 +76,7 @@ double Recon::estimateOffset(const MeasurementInfo aExpInfo, const CEInfo& aCe,
// Digitalfilterdelay=-(4.67e-6 + 0.8e-6);
// %seconds
// elseif(strcmp(info.Hardware,'USCT3Dv3')==1)
double digitalfilterdelay = 0;
float digitalfilterdelay = 0;
// else
// Digitalfilterdelay = 0;
// end

View File

@@ -5,7 +5,7 @@
namespace Recon
{
double estimateOffset(const MeasurementInfo aExpInfo, const CEInfo& aCe, const Aurora::Matrix& aMatchedFilter);
float estimateOffset(const MeasurementInfo aExpInfo, const CEInfo& aCe, const Aurora::Matrix& aMatchedFilter);
}
#endif

View File

@@ -12,11 +12,11 @@ Matrix Recon::imageExtrapolation(const Matrix& aImg, int aDim)
if(aDim == 3)
{
size_t zSize = aImg.getDimSize(2);
Matrix zIndexZero = Matrix::fromRawData(new double[zSize], zSize);
Matrix zIndexZero = Matrix::fromRawData(new float[zSize], zSize);
std::vector<int> indexe;
for(int i=0; i<zSize; ++i)
{
double value = sum(aImg($,$,i).toMatrix(), FunctionDirection::All)[0];
float value = sum(aImg($,$,i).toMatrix(), FunctionDirection::All)[0];
if(value == 0)
{
zIndexZero[i] = 1;

View File

@@ -18,7 +18,7 @@ void Recon::precalcImageParameters(const GeometryInfo& aGeom)
minSize.forceReshape(2, 1, 1);
reflectParams::imageResolution = max(maxSize - minSize, FunctionDirection::All)[0] / reflectParams::pixelResolutionX;
reflectParams::pixelResolutionZ = std::floor((aGeom.maxSize[2] - aGeom.minSize[2]) / reflectParams::imageResolution) + 1;
double* imageXYZData = Aurora::malloc(3);
float* imageXYZData = Aurora::malloc(3);
imageXYZData[0] = ceil(reflectParams::pixelResolutionX * (reflectParams::imageEndpoint[0] - reflectParams::imageStartpoint[0]) /
(aGeom.maxSize[0] - aGeom.minSize[0]));
imageXYZData[1] = ceil(reflectParams::pixelResolutionY * (reflectParams::imageEndpoint[1] - reflectParams::imageStartpoint[1]) /

View File

@@ -18,7 +18,7 @@
#include <vector>
namespace Recon {
Aurora::Matrix _createDiffVector(const Aurora::Matrix & aMatrix, double beginValue,double endValue){
Aurora::Matrix _createDiffVector(const Aurora::Matrix & aMatrix, float beginValue,float endValue){
auto r = Aurora::zeros(1,aMatrix.getDataSize()+1);
r[0] = beginValue;
for (size_t i = 1; i < aMatrix.getDataSize(); i++)
@@ -46,8 +46,8 @@ namespace Recon {
{
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);
if (op == 1)vsSubI(result.getDataSize(),aMatrix.getData() , 1, aMatrix.getData()+1,1,result.getData(),1);
if (op == 0)vsAddI(result.getDataSize(),aMatrix.getData()+1 , 1, aMatrix.getData(),1,result.getData(),1);
return result;
}
else{
@@ -77,7 +77,7 @@ 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};
std::shared_ptr<float> ptr{aMatrix.getData()+aMatrix.getDimSize(0)*col*aMatrix.getValueType(),fakeFree};
return Aurora::Matrix(ptr,{aMatrix.getDimSize(0),1,1},aMatrix.getValueType());
}
@@ -93,8 +93,8 @@ namespace Recon {
Aurora::Matrix offsetSignalFourier(const PreComputes &preComputes, int nfft)
{
double dt = -(preComputes.offset/preComputes.timeInterval);
double binStart = floor((double)nfft/2);
float dt = -(preComputes.offset/preComputes.timeInterval);
float binStart = floor((float)nfft/2);
auto temp = Aurora::zeros(1,nfft);
for (size_t i = 0; i < temp.getDataSize(); i++)
{
@@ -105,7 +105,7 @@ namespace Recon {
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);
cblas_scopy(fftBin.getDataSize(), fftBin.getData(), 1, negCfftbin.getData()+1, 2);
auto exponent = exp(negCfftbin);
if (!exponent.isComplex())exponent = complex(exponent);
return exponent;
@@ -170,7 +170,7 @@ namespace Recon {
if (blockedSL[i] == blockedRL[i])
{
auto begin_ptr = _blockedAScans.getData();
int count = reflectParams::suppressSameHeadLength* sizeof(double);
int count = reflectParams::suppressSameHeadLength* sizeof(float);
std::memset(begin_ptr,0,count);
}
}
@@ -192,23 +192,23 @@ namespace Recon {
else{
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);
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);
float* value1 = Aurora::malloc(fx.getDataSize());
vsMulI(fx.getDataSize(), fx.getData(), 2, fh.getData(), 2, value1, 1);
float* value2 = Aurora::malloc(fx.getDataSize());
vsMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData() + 1, 2, value2, 1);
float* realData = Aurora::malloc(fx.getDataSize());
vsAdd(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);
vsMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData(), 2, value1, 1);
vsMulI(fx.getDataSize(), fx.getData(), 2, fh.getData() + 1, 2, value2, 1);
float* imagData = Aurora::malloc(fx.getDataSize());
vsSub(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);
float* complexData = Aurora::malloc(real.getDataSize(), true);
cblas_scopy(real.getDataSize(), real.getData(), 1 , complexData ,2);
cblas_scopy(image.getDataSize(), image.getData(), 1 , complexData + 1 ,2);
Aurora::Matrix complex = Aurora::Matrix::New(complexData, real.getDimSize(0), real.getDimSize(1), 1, Aurora::Complex);
Aurora::free(value1);
Aurora::free(value2);
@@ -267,9 +267,9 @@ namespace Recon {
auto mean_energySum = Aurora::zeros(l_partData + windowLength,1);
for (size_t k = 0; k < windowLength; k++)
{
double* begin = mean_energySum.getData() + k;
float* begin = mean_energySum.getData() + k;
int length = l_partData;
vdAddI(length, energySum.getData(), 1, begin, 1, begin, 1);
vsAddI(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));
@@ -288,7 +288,7 @@ namespace Recon {
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();
float median_mean_energySum = Aurora::median(mean_energySum).getScalar();
bool flag = false;
if (mean_energySum[index] > median_mean_energySum)
{
@@ -319,16 +319,16 @@ namespace Recon {
int winLength2 = 100;
winLength2 = sig_end+winLength2+200-sig_begin+1;
auto win = generateGaussWindow(winLength2);
double* begin = _blockedAScans.getData() + (size_t)sig_begin-1;
float* begin = _blockedAScans.getData() + (size_t)sig_begin-1;
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);
vsMulI(length, begin, 1, winReverse.getData(), 1, begin, 1);
winLength2 = sig_end;
win = generateGaussWindow(winLength2);
begin = _blockedAScans.getData() ;
length = sig_end;
vdMulI(length, begin, 1, win.getData(), 1, begin, 1);
vsMulI(length, begin, 1, win.getData(), 1, begin, 1);
// _blockedAScans = fft(_blockedAScans);
// blockedAScans(Aurora::$,i) =_blockedAScans;
}
@@ -338,14 +338,14 @@ namespace Recon {
blockedAScans = fft(blockedAScans);
if(reflectParams::useOptPulse==1){
auto n = nSamples;
auto nHalf = round((double)n/2);
auto nHalf = round((float)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(0,Aurora::$) = std::complex<float>{0,0};
_blockedAScans.setBlockComplexValue(0,nHalf,_blockedAScans.getDimSize(0)-1, std::complex<float>{0,0});
_blockedAScans.setBlock(0,1,nHalf-1,_blockedAScans.block(0,1,nHalf-1)*2);
}
@@ -487,9 +487,9 @@ namespace Recon {
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();
float highNoiseScoreMean = mean(highNoiseScore).getScalar();
float highNoiseScoreStd = Aurora::std(highNoiseScore).getScalar();
float treshold99= (mean(highNoiseScore)+2.596*Aurora::std(highNoiseScore)).getScalar();
Aurora::compareSet(valid,ascanMapValue,treshold99,0,Aurora::GT);
// debugOutput.blockStatisticValue = ascanMapValue;
// debugOutput.blockStatisticUsedData = ascanMapValue > treshold99;

View File

@@ -24,9 +24,9 @@ PreprocessReflectionData Recon::preprocessTransmissionReconstructionForReflectio
Matrix maxPos;
if(!aDdims.isNull())
{
minPos = Matrix::fromRawData(new double[3], 1, 3);
minPos = Matrix::fromRawData(new float[3], 1, 3);
size_t maxPosSize = aDdims.getDataSize() - 3;
maxPos = Matrix::fromRawData(new double[maxPosSize], 1, maxPosSize);
maxPos = Matrix::fromRawData(new float[maxPosSize], 1, maxPosSize);
std::copy(aDdims.getData(), aDdims.getData() + 3, minPos.getData());
std::copy(aDdims.getData() + 3, aDdims.getData() + 3 + maxPosSize, maxPos.getData());
}
@@ -96,7 +96,7 @@ PreprocessReflectionData Recon::preprocessTransmissionReconstructionForReflectio
if(soundSpeedCorrection == 0 && attenuationCorrection == 0)
{
result.saftMode = 2;
double deltaTransMap = reflectParams::resolutionTransmMap;
float deltaTransMap = reflectParams::resolutionTransmMap;
Matrix minTransdPos = min(aGeometryInfo.minEmitter, aGeometryInfo.minReceiver);
Matrix maxTransdPos = max(aGeometryInfo.maxEmitter, aGeometryInfo.maxReceiver);

View File

@@ -12,15 +12,15 @@ namespace Recon
Aurora::Matrix speedMap3D;
Aurora::Matrix attenuationMap3D;
Aurora::Matrix beginTransMap;
double deltaTransMap;
float deltaTransMap;
};
struct PreprocessReflectionData
{
TransRecos transRecos;
int saftMode;
double attenuationCorrection;
double soundSpeedCorrection;
float attenuationCorrection;
float soundSpeedCorrection;
};
PreprocessReflectionData preprocessTransmissionReconstructionForReflection(const Aurora::Matrix& aSOSMap, const Aurora::Matrix& aATTMap,
const Aurora::Matrix& aDdims, const GeometryInfo& aGeometryInfo,

View File

@@ -17,9 +17,9 @@ TransmissionVolme Recon::resampleTransmissionVolume(const Aurora::Matrix& aTrans
const Aurora::Matrix& aMinPos,const GeometryInfo& aGeom)
{
TransmissionVolme result;
double resSOSX = (aMaxPos[0] - aMinPos[0]) / (aTransVol.getDimSize(0) - 1);
double resSOSY = (aMaxPos[1] - aMinPos[1]) / (aTransVol.getDimSize(1) - 1);
double resSOSZ = (aMaxPos[2] - aMinPos[2]) / (aTransVol.getDimSize(2) - 1);
float resSOSX = (aMaxPos[0] - aMinPos[0]) / (aTransVol.getDimSize(0) - 1);
float resSOSY = (aMaxPos[1] - aMinPos[1]) / (aTransVol.getDimSize(1) - 1);
float resSOSZ = (aMaxPos[2] - aMinPos[2]) / (aTransVol.getDimSize(2) - 1);
Matrix xSOS = createVectorMatrix(aMinPos[0] + 0.5 * resSOSX, resSOSX, aMaxPos[0] + 0.5 * resSOSX);
Matrix ySOS = createVectorMatrix(aMinPos[1] + 0.5 * resSOSY, resSOSY, aMaxPos[1] + 0.5 * resSOSY);
@@ -30,7 +30,7 @@ TransmissionVolme Recon::resampleTransmissionVolume(const Aurora::Matrix& aTrans
Matrix minTransdPos = min(aGeom.minEmitter, aGeom.minReceiver);
Matrix maxTransdPos = max(aGeom.maxEmitter, aGeom.maxReceiver);
double deltaTransMap = (aMaxPos[0] - aMinPos[0]) / (aTransVol.getDimSize(0) - 1);
float deltaTransMap = (aMaxPos[0] - aMinPos[0]) / (aTransVol.getDimSize(0) - 1);
Matrix beginTransMap = transpose(min(transpose(reflectParams::imageStartpoint), minTransdPos)) - 3 * reflectParams::resolutionTransmMap;
Matrix endSpeedMap = max(transpose(reflectParams::imageEndpoint), maxTransdPos) + 3 * reflectParams::resolutionTransmMap;
endSpeedMap[2] = endSpeedMap[2] + 0.025;

View File

@@ -9,7 +9,7 @@ namespace Recon
struct TransmissionVolme
{
Aurora::Matrix transMap;
double deltaTransMap;
float deltaTransMap;
Aurora::Matrix beginTransMap;
};

View File

@@ -42,7 +42,7 @@ namespace Recon {
const Aurora::Matrix &senderIndex,
const Aurora::Matrix &receiverPosGeom,
const Aurora::Matrix &senderPosGeom,
int SAFT_mode, double TimeInterval,
int SAFT_mode, float TimeInterval,
const TransRecos& transRecos,
const Aurora::Matrix &Env)
{
@@ -228,7 +228,7 @@ namespace Recon {
// delete [] fdata4;
// delete [] fdata5;
// delete [] gdata;
// return Aurora::Matrix::fromRawData((double *)result.Data, result.Dims[0],
// return Aurora::Matrix::fromRawData((float *)result.Data, result.Dims[0],
// result.Dims[1], result.Dims[2]);
return Aurora::Matrix();
@@ -245,7 +245,7 @@ namespace Recon {
delete [] fdata4;
delete [] fdata5;
delete [] gdata;
return Aurora::Matrix::fromRawData((double *)result.Data, result.Dims[0],
return Aurora::Matrix::fromRawData((float *)result.Data, result.Dims[0],
result.Dims[1], result.Dims[2]);
}
}
@@ -254,7 +254,7 @@ namespace Recon {
recontructSAFT(const Aurora::Matrix &AScans, const Aurora::Matrix &senderPos,
const Aurora::Matrix &receiverPos, const Aurora::Matrix &mpIndex,
int SAFT_mode, TransRecos &transRecos, Aurora::Matrix &Env) {
double TimeInterval = 1.0 / (reflectParams::aScanReconstructionFrequency);
float TimeInterval = 1.0 / (reflectParams::aScanReconstructionFrequency);
std::vector<int> motorPosAvailable;
std::unique_copy(mpIndex.getData(), mpIndex.getData() + mpIndex.getDataSize(),
std::back_inserter(motorPosAvailable));

View File

@@ -69,13 +69,13 @@ Aurora::Matrix Recon::startReflectionReconstruction( Parser* aParser, int aSAFT_
Matrix sn = aSnList.block(0, transParams::senderElementSize*k, transParams::senderElementSize*k+transParams::senderElementSize - 1);
auto blockData = getAscanBlockPreprocessed(aParser, mp, sl, sn, aRlList, aRnList, aGeom, aExpInfo, true, false);
double* channelListSizeData = Aurora::malloc(2);
float* channelListSizeData = Aurora::malloc(2);
channelListSizeData[0] = channelList.getDimSize(0);
channelListSizeData[1] = channelList.getDimSize(1);
Matrix channelListSize = Matrix::New(channelListSizeData, 2, 1);
Matrix ind = sub2ind(channelListSize, {blockData.rlBlock, blockData.rnBlock});
size_t channelBlockSize = ind.getDataSize();
double* channelBlockData = Aurora::malloc(channelBlockSize);
float* channelBlockData = Aurora::malloc(channelBlockSize);
for(size_t i=0; i<channelBlockSize; ++i)
{
channelBlockData[i] = channelList[ind[i] - 1];