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

commit dev-sun.

Browse Source
This commit is contained in:
sunwen
2023-12-22 11:17:18 +08:00
parent d015b38845
commit 410d657fe7
27 changed files with 1346 additions and 267 deletions
Download Patch File Download Diff File Expand all files Collapse all files

253
src/transmissionReconstruction/detection/getTransmissionData.cpp
View File

@@ -1,4 +1,7 @@
#include "getTransmissionData.h"
#include "getTransmissionData.cuh"
#include "AuroraDefs.h"
#include "CudaMatrix.h"
#include "Function.h"
#include "Function1D.h"
#include "Function2D.h"
@@ -28,6 +31,11 @@
#include <mutex>
#include <condition_variable>
#include <cuda_runtime.h>
#include "Function1D.cuh"
#include "Function2D.cuh"
#include <sys/time.h>
using namespace Recon;
using namespace Aurora;
@@ -42,8 +50,8 @@ namespace
Matrix waterTempBlock;
MetaInfos metaInfos;
Matrix ascanBlock;
Matrix ascanBlockRef;
CudaMatrix ascanBlock;
CudaMatrix ascanBlockRef;
Matrix dists;
Matrix distRefBlock;
Matrix waterTempRefBlock;
@@ -56,28 +64,59 @@ namespace
std::mutex PROCESS_BUFFER_MUTEX;
std::condition_variable PROCESS_BUFFER_CONDITION;
int BUFFER_COUNT = 0;
int BUFFER_SIZE = 3;
int BUFFER_SIZE = 4;//<=8
Matrix prepareAScansForTransmissionDetection(const Matrix& aAscanBlock, const Matrix& aGainBlock)
void printTime()
{
struct timeval tpend;
gettimeofday(&tpend,NULL);
int secofday = (tpend.tv_sec + 3600 * 8 ) % 86400;
int hours = secofday / 3600;
int minutes = (secofday - hours * 3600 ) / 60;
int seconds = secofday % 60;
int milliseconds = tpend.tv_usec/1000;
std::cout<< hours << ":" <<minutes<<":"<<seconds<<"."<<milliseconds<<std::endl;
}
CudaMatrix prepareAScansForTransmissionDetection(const CudaMatrix& aAscanBlock, const CudaMatrix& aGainBlock)
{
Matrix result = aAscanBlock / repmat(aGainBlock, aAscanBlock.getDimSize(0), 1);
CudaMatrix result = aAscanBlock / repmat(aGainBlock, aAscanBlock.getDimSize(0), 1);
result = result - repmat(mean(result,FunctionDirection::Column), result.getDimSize(0), 1);
return result;
}
Aurora::CudaMatrix calculateSnr(const Aurora::CudaMatrix &aMDataBlock,
float aReferenceNoise) {
auto maxSignal = max(abs(aMDataBlock));
auto snrBlock = 10 * log(maxSignal / aReferenceNoise, 10);
return snrBlock;
}
BlockOfTransmissionData getBlockOfTransmissionData(const Matrix& aMp, const Matrix& aMpRef, const Matrix& aSl, const Matrix& aSn, const Matrix& aRlList, const Matrix& aRnList,
const TasTemps& aTasTemps, const Matrix& aExpectedSOSWater, GeometryInfo aGeom, GeometryInfo& aGeomRef,
const Matrix& aSnrRmsNoise, const Matrix& aSnrRmsNoiseRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
{
BlockOfTransmissionData result;
MetaInfos metaInfos;
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);
Matrix ascanBlock = prepareAScansForTransmissionDetection(blockData.ascanBlockPreprocessed, blockData.gainBlock);
Matrix ascanBlockRef = prepareAScansForTransmissionDetection(blockDataRef.ascanBlockPreprocessed, blockDataRef.gainBlock);
blockData.ascanBlockPreprocessed = Matrix();
blockDataRef.ascanBlockPreprocessed = Matrix();
MetaInfos metaInfos;
//printTime();
//2500ms
auto blockData = getAscanBlockPreprocessedCuda(aParser, aMp, aSl, aSn, aRlList, aRnList, aGeom, aExpInfo, true, true);
auto blockDataRef = getAscanBlockPreprocessedCuda(aParserRef, aMpRef, aSl, aSn, aRlList, aRnList, aGeomRef, aExpInfoRef, true, true);
//printTime();
//180ms
// auto t1 = blockData.ascanBlockPreprocessed.toDeviceMatrix();
// auto t2 = blockDataRef.ascanBlockPreprocessed.toDeviceMatrix();
// auto t3 = blockData.gainBlock.toDeviceMatrix();
// auto t4 = blockDataRef.gainBlock.toDeviceMatrix();
//printTime();
//20ms
CudaMatrix ascanBlock = prepareAScansForTransmissionDetection(blockData.ascanBlockPreprocessed,blockData.gainBlock);
CudaMatrix ascanBlockRef = prepareAScansForTransmissionDetection(blockDataRef.ascanBlockPreprocessed,blockDataRef.gainBlock);
//printTime();
//20ms
blockData.ascanBlockPreprocessed = CudaMatrix();
blockDataRef.ascanBlockPreprocessed = CudaMatrix();
if(aExpInfo.Hardware == "USCT3dv3")
{
Matrix channelList = precalculateChannelList(aRlList, aRnList, aExpInfo, aPreComputes);
@@ -93,71 +132,48 @@ namespace
channelListBlockData[i] = channelList[ind[i] - 1];
}
Matrix channelListBlock = Matrix::New(channelListBlockData, 1, channelListBlockSize);
Matrix fx = fft(ascanBlock);
float* fhData = Aurora::malloc(aExpInfo.matchedFilter.getDimSize(0) * channelListBlockSize, true);
Matrix fh = Matrix::New(fhData, aExpInfo.matchedFilter.getDimSize(0), channelListBlockSize, 1, Aurora::Complex);
size_t matchedFilterRowDataSize = aExpInfo.matchedFilter.getDimSize(0)*2;
for(size_t i=0; i<channelListBlockSize; ++i)
{
cblas_scopy(matchedFilterRowDataSize, aExpInfo.matchedFilter.getData() + (size_t)(channelListBlock[i] - 1) * matchedFilterRowDataSize, 1 , fhData ,1);
fhData += matchedFilterRowDataSize;
}
// Matrix fxReal = Aurora::real(fx);
// Matrix fhReal = Aurora::real(fh);
// Matrix fxImag = Aurora::imag(fx);
// Matrix fhImag = Aurora::imag(fh);
// Matrix real = fxReal * fhReal + fxImag * fhImag;
// Matrix image = fxImag * fhReal - fxReal * fhImag;
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);
Matrix real = Matrix::New(realData, fx.getDimSize(0), fx.getDimSize(1));
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);
Matrix image = Matrix::New(imagData, fx.getDimSize(0), fx.getDimSize(1));
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);
Matrix complex = Matrix::New(complexData, real.getDimSize(0), real.getDimSize(1), 1, Aurora::Complex);
//printTime();
//20ms
CudaMatrix fx = fft(ascanBlock);
//printTime();
//50ms
// float* fhData = nullptr;
// cudaMalloc((void**)&fhData, sizeof(float) * aExpInfo.matchedFilter.getDimSize(0) * channelListBlockSize * Aurora::Complex);
// CudaMatrix fh = CudaMatrix::fromRawData(fhData, aExpInfo.matchedFilter.getDimSize(0), channelListBlockSize, 1, Aurora::Complex);
// size_t matchedFilterRowDataSize = aExpInfo.matchedFilter.getDimSize(0)*2;
// for(size_t i=0; i<channelListBlockSize; ++i)
// {
// cudaMemcpy(fhData, aExpInfo.matchedFilter.getData() + (size_t)(channelListBlock[i] - 1) * matchedFilterRowDataSize, sizeof(float) * matchedFilterRowDataSize, cudaMemcpyHostToDevice);
// fhData += matchedFilterRowDataSize;
// }
Aurora::CudaMatrix matchedFilterDevice = aExpInfo.matchedFilter.toDeviceMatrix();
Aurora::CudaMatrix channelListBlockDevice = channelListBlock.toDeviceMatrix();
Aurora::CudaMatrix fh = createFhMatrix(matchedFilterDevice, channelListBlockDevice);
//printTime();
//20ms
CudaMatrix complex = getTransmissionDataSubFunction(fx, fh);
ascanBlock = Aurora::real(ifft(complex));
//printTime();
//20s
fx = fft(ascanBlockRef);
fhData = Aurora::malloc(aExpInfoRef.matchedFilter.getDimSize(0) * channelListBlockSize, true);
fh = Matrix::New(fhData, aExpInfoRef.matchedFilter.getDimSize(0), channelListBlockSize, 1, Aurora::Complex);
matchedFilterRowDataSize = aExpInfoRef.matchedFilter.getDimSize(0)*2;
for(size_t i=0; i<channelListBlockSize; ++i)
{
cblas_scopy(matchedFilterRowDataSize, aExpInfoRef.matchedFilter.getData() + (size_t)(channelListBlock[i] - 1) * matchedFilterRowDataSize, 1 , fhData ,1);
fhData += matchedFilterRowDataSize;
}
// real = Aurora::real(fx) * Aurora::real(fh) + Aurora::imag(fx) * Aurora::imag(fh);
// image = Aurora::imag(fx) * Aurora::real(fh) - Aurora::real(fx) * Aurora::imag(fh);
vsMulI(fx.getDataSize(), fx.getData(), 2, fh.getData(), 2, value1, 1);
vsMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData() + 1, 2, value2, 1);
realData = Aurora::malloc(fx.getDataSize());
vsAdd(fx.getDataSize(), value1, value2, realData);
real = Matrix::New(realData, fx.getDimSize(0), fx.getDimSize(1));
vsMulI(fx.getDataSize(), fx.getData() + 1, 2, fh.getData(), 2, value1, 1);
vsMulI(fx.getDataSize(), fx.getData(), 2, fh.getData() + 1, 2, value2, 1);
imagData = Aurora::malloc(fx.getDataSize());
vsSub(fx.getDataSize(), value1, value2, imagData);
image = Matrix::New(imagData, fx.getDimSize(0), fx.getDimSize(1));
Aurora::free(value1);
Aurora::free(value2);
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);
complex = Matrix::New(complexData, real.getDimSize(0), real.getDimSize(1), 1, Aurora::Complex);
//printTime();
//50ms
// cudaMalloc((void**)&fhData, sizeof(float) * aExpInfo.matchedFilter.getDimSize(0) * channelListBlockSize * Aurora::Complex);
// fh = CudaMatrix::fromRawData(fhData, aExpInfoRef.matchedFilter.getDimSize(0), channelListBlockSize, 1, Aurora::Complex);
// matchedFilterRowDataSize = aExpInfoRef.matchedFilter.getDimSize(0)*2;
// for(size_t i=0; i<channelListBlockSize; ++i)
// {
// cudaMemcpy(fhData, aExpInfoRef.matchedFilter.getData() + (size_t)(channelListBlock[i] - 1) * matchedFilterRowDataSize, sizeof(float) * matchedFilterRowDataSize, cudaMemcpyHostToDevice);
// fhData += matchedFilterRowDataSize;
// }
matchedFilterDevice = aExpInfoRef.matchedFilter.toDeviceMatrix();
fh = createFhMatrix(matchedFilterDevice, channelListBlockDevice);
//printTime();
//20ms
complex = getTransmissionDataSubFunction(fx, fh);
ascanBlockRef = Aurora::real(ifft(complex));
//printTime();
//20ms
}
else
{
@@ -167,24 +183,27 @@ namespace
if(transParams::applyCalib)
{
metaInfos.snrValues = calculateSnr(ascanBlock, aSnrRmsNoise[0]);
metaInfos.snrValuesRef = calculateSnr(ascanBlockRef, aSnrRmsNoiseRef[0]);
metaInfos.snrValues = calculateSnr(ascanBlock, aSnrRmsNoise[0]).toHostMatrix();
metaInfos.snrValuesRef = calculateSnr(ascanBlockRef, aSnrRmsNoiseRef[0]).toHostMatrix();
}
// printTime();
//3ms
Matrix dists = distanceBetweenTwoPoints(blockData.senderPositionBlock, blockData.receiverPositionBlock);
Matrix distRefBlock = distanceBetweenTwoPoints(blockDataRef.senderPositionBlock, blockDataRef.receiverPositionBlock);
//printTime();
//2ms
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;
metaInfos.slBlock = blockData.slBlock;
metaInfos.snBlock = blockData.snBlock;
metaInfos.rlBlock = blockData.rlBlock;
metaInfos.rnBlock = blockData.rnBlock;
}
// printTime();
// 1ms
// if(transParams::saveDetection || transParams::outlierOnTasDetection || transParams::saveDebugInfomation)
// {
// metaInfos.mpBlock = blockData.mpBlock;
// metaInfos.slBlock = blockData.slBlock;
// metaInfos.snBlock = blockData.snBlock;
// metaInfos.rlBlock = blockData.rlBlock;
// metaInfos.rnBlock = blockData.rnBlock;
// }
result.metaInfos = metaInfos;
result.senderBlock = blockData.senderPositionBlock;
result.receiverBlock = blockData.receiverPositionBlock;
@@ -199,7 +218,7 @@ namespace
// DetectResult detect = transmissionDetection(ascanBlock, ascanBlockRef, dists, distRefBlock, waterTempBlock, waterTempRefBlock, aExpectedSOSWater[0]);
// result.attData = detect.att;
// result.tofData = detect.tof;
//printTime();
return result;
}
@@ -208,8 +227,9 @@ namespace
void getBlockOfTransmissionDataInThread(size_t aIndex, const Matrix& aMp, const Matrix& aMpRef, const Matrix& aSl, const Matrix& aSn, const Matrix& aRlList, const Matrix& aRnList,
const TasTemps& aTasTemps, const Matrix& aExpectedSOSWater, GeometryInfo aGeom, GeometryInfo aGeomRef,
const Matrix& aSnrRmsNoise, const Matrix& aSnrRmsNoiseRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef, unsigned int aGPUId)
{
cudaSetDevice(aGPUId);
auto buffer = getBlockOfTransmissionData(aMp, aMpRef, aSl, aSn, aRlList, aRnList, aTasTemps,
aExpectedSOSWater, aGeom, aGeomRef, aSnrRmsNoise, aSnrRmsNoiseRef,
aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef);
@@ -245,7 +265,7 @@ void createThreadForGetBlockOfTransmissionData(const Matrix& aMotorPos, const M
CREATE_BUFFER_CONDITION.wait(lock, []{return BUFFER_COUNT<BUFFER_SIZE;});
++BUFFER_COUNT;
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);
speedUpThread[index] = std::thread(getBlockOfTransmissionDataInThread,index,mp,mpRef,sl,sn,aRlList,aRnList,aTasTemps,aExpectedSOSWater,aGeom,aGeomRef,aSnrRmsNoise,aSnrRmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef, index % BUFFER_SIZE);
}
}
}
@@ -334,8 +354,9 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
lock.unlock();
auto blockData = BLOCK_OF_TRANSIMISSIONDARA_BUFFER[std::to_string(index)];
cudaSetDevice(index % BUFFER_SIZE);
DetectResult detect = transmissionDetection( blockData.ascanBlock, blockData.ascanBlockRef,
blockData.dists, blockData.distRefBlock,
blockData.dists.toDeviceMatrix(), blockData.distRefBlock.toDeviceMatrix(),
blockData.waterTempBlock, blockData.waterTempRefBlock,
aTemp.expectedSOSWater[0]);
blockData.attData = detect.att;
@@ -355,14 +376,14 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
cblas_scopy(numUsedData, transmissionBlock.attData.getData(), 1, attDataTotal.getData() + numData, 1);
cblas_scopy(numUsedData, transmissionBlock.waterTempBlock.getData(), 1, waterTempList.getData() + numData, 1);
if(transParams::saveDetection || transParams::outlierOnTasDetection || transParams::saveDebugInfomation)
{
cblas_scopy(numUsedData, transmissionBlock.metaInfos.mpBlock.getData(), 1, mpBlockTotal.getData() + numData, 1);
cblas_scopy(numUsedData, transmissionBlock.metaInfos.slBlock.getData(), 1, slBlockTotal.getData() + numData, 1);
cblas_scopy(numUsedData, transmissionBlock.metaInfos.snBlock.getData(), 1, snBlockTotal.getData() + numData, 1);
cblas_scopy(numUsedData, transmissionBlock.metaInfos.rlBlock.getData(), 1, rlBlockTotal.getData() + numData, 1);
cblas_scopy(numUsedData, transmissionBlock.metaInfos.rnBlock.getData(), 1, rnBlockTotal.getData() + numData, 1);
}
// if(transParams::saveDetection || transParams::outlierOnTasDetection || transParams::saveDebugInfomation)
// {
// cblas_scopy(numUsedData, transmissionBlock.metaInfos.mpBlock.getData(), 1, mpBlockTotal.getData() + numData, 1);
// cblas_scopy(numUsedData, transmissionBlock.metaInfos.slBlock.getData(), 1, slBlockTotal.getData() + numData, 1);
// cblas_scopy(numUsedData, transmissionBlock.metaInfos.snBlock.getData(), 1, snBlockTotal.getData() + numData, 1);
// cblas_scopy(numUsedData, transmissionBlock.metaInfos.rlBlock.getData(), 1, rlBlockTotal.getData() + numData, 1);
// cblas_scopy(numUsedData, transmissionBlock.metaInfos.rnBlock.getData(), 1, rnBlockTotal.getData() + numData, 1);
// }
numData += numUsedData;
std::unique_lock<std::mutex> lockBufferCount(CREATE_BUFFER_MUTEX);
BLOCK_OF_TRANSIMISSIONDARA_BUFFER.erase(std::to_string(index));
@@ -371,7 +392,7 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
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))));
//Recon::notifyProgress(6+10*((i+1)*(j+1)/(aMotorPos.getDataSize()*(aSlList.getDataSize()/ transParams::senderTASSize))));
}
}
speedUpThread.join();
@@ -399,14 +420,14 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
tofDataTotal = removeDataFromArrays(tofDataTotal, filter);
attDataTotal = removeDataFromArrays(attDataTotal, filter);
waterTempList = removeDataFromArrays(waterTempList, filter);
if(transParams::saveDebugInfomation || transParams::outlierOnTasDetection || transParams::saveDetection)
{
mpBlockTotal = removeDataFromArrays(mpBlockTotal, filter);
slBlockTotal = removeDataFromArrays(slBlockTotal, filter);
snBlockTotal = removeDataFromArrays(snBlockTotal, filter);
rlBlockTotal = removeDataFromArrays(rlBlockTotal, filter);
rnBlockTotal = removeDataFromArrays(rnBlockTotal, filter);
}
// if(transParams::saveDebugInfomation || transParams::outlierOnTasDetection || transParams::saveDetection)
// {
// mpBlockTotal = removeDataFromArrays(mpBlockTotal, filter);
// slBlockTotal = removeDataFromArrays(slBlockTotal, filter);
// snBlockTotal = removeDataFromArrays(snBlockTotal, filter);
// rlBlockTotal = removeDataFromArrays(rlBlockTotal, filter);
// rnBlockTotal = removeDataFromArrays(rnBlockTotal, filter);
// }
Matrix valid;
if(transParams::applyCalib)
@@ -432,14 +453,14 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
}
dataInfno.findDefect = Matrix::New(findDefectData, 1, findDefectDataIndex);
if(transParams::saveDebugInfomation)
{
dataInfno.sn = snBlockTotal;
dataInfno.sl = slBlockTotal;
dataInfno.rn = rnBlockTotal;
dataInfno.rl = rlBlockTotal;
dataInfno.mp = mpBlockTotal;
}
// if(transParams::saveDebugInfomation)
// {
// dataInfno.sn = snBlockTotal;
// dataInfno.sl = slBlockTotal;
// dataInfno.rn = rnBlockTotal;
// dataInfno.rl = rlBlockTotal;
// dataInfno.mp = mpBlockTotal;
// }
tofDataTotal = removeDataFromArrays(tofDataTotal, valid);
attDataTotal = removeDataFromArrays(attDataTotal, valid);