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

Optimize transmission preprocess and adapt new cuda version.

Browse Source
This commit is contained in:
sunwen
2024-12-19 10:07:47 +08:00
parent af865bb815
commit 69d95e8de9
19 changed files with 147 additions and 602 deletions
Download Patch File Download Diff File Expand all files Collapse all files

92
src/Aurora.cu
View File

@@ -66,23 +66,6 @@ __global__ void validKernel(const float* aData, const float* aValid, float* aOut
}
}
// __global__ void validSubKernel(const double* aValid, double* aOutput, unsigned int* aCount, int aValidSize)
// {
// int index = blockIdx.x * blockDim.x + threadIdx.x;
// if(index == 0)
// {
// for(int i=0;i<aValidSize;++i)
// {
// if(aValid[i] == 1)
// {
// aOutput[*aCount] = i;
// ++(*aCount);
// }
// }
// }
// __syncthreads();
// }
Aurora::CudaMatrix Aurora::valid(const Aurora::CudaMatrix aData, const Aurora::CudaMatrix aValid)
{
int validSize = aValid.getDataSize();
@@ -116,81 +99,6 @@ Aurora::CudaMatrix Aurora::valid(const Aurora::CudaMatrix aData, const Aurora::C
return Aurora::CudaMatrix::fromRawData(result, rowCount, validColumnCount);
}
texture<float, cudaTextureType2D, cudaReadModeElementType> tex;
cudaArray* array;
__global__ void testKernel(float* aData,cudaTextureObject_t aTexObj, cudaSurfaceObject_t aSurface)
{
float a = tex2D(tex,5.5,5.5);
float b = tex2D<float>(aTexObj,5.5,5.5);
float2 c = tex2D<float2>(aSurface,1,1);
printf("%f\n",a);
printf("%f\n",b);
printf("%f\n",c.x);
printf("%f\n",c.y);
}
__global__ void writeSurfaceKernel( cudaSurfaceObject_t aSurface)
{
float2 value;
value.x = 100;
value.y = 99;
surf2Dwrite(value, aSurface, 1, 1 );
}
void subTest(cudaTextureObject_t& aTexture)
{
cudaResourceDesc resourceDesc;
cudaTextureDesc textureDesc;
memset(&resourceDesc, 0, sizeof(resourceDesc));
resourceDesc.resType = cudaResourceTypeArray;
resourceDesc.res.array.array = array; // 指向设备端的 CUDA 数组
// 在 textureDesc 中设置纹理描述
memset(&textureDesc, 0, sizeof(textureDesc));
textureDesc.addressMode[0] = cudaAddressModeClamp;
textureDesc.addressMode[1] = cudaAddressModeClamp;
textureDesc.filterMode = cudaFilterModeLinear;
textureDesc.readMode = cudaReadModeElementType;
textureDesc.normalizedCoords = false;
//textureDesc.channelDesc = texChannelDescSpeedOfSoundField;
cudaCreateTextureObject(&aTexture, &resourceDesc, &textureDesc, nullptr);
}
void Aurora::test(float* aData)
{
tex.addressMode[0] = cudaAddressModeClamp; // Texturreferenz beschreiben
tex.addressMode[1] = cudaAddressModeClamp;
tex.filterMode = cudaFilterModeLinear;
tex.normalized = 0;
cudaChannelFormatDesc texChannelDescSpeedOfSoundField = cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindFloat);
cudaMallocArray(&array, &texChannelDescSpeedOfSoundField, 10, 9);
cudaMemcpyToArray(array, 0, 0, aData,10 * 9 *sizeof(float), cudaMemcpyHostToDevice);
cudaBindTextureToArray ( &tex, array, &texChannelDescSpeedOfSoundField );
cudaTextureObject_t textureObj;
subTest(textureObj);
struct cudaResourceDesc resDesc;
memset(&resDesc, 0, sizeof(resDesc));
resDesc.resType = cudaResourceTypeArray;
// Create the surface objects
resDesc.res.array.array = array;
cudaSurfaceObject_t inputSurfObj = 0;
cudaCreateSurfaceObject(&inputSurfObj, &resDesc);
writeSurfaceKernel<<<1,1>>>(inputSurfObj);
cudaDeviceSynchronize();
testKernel<<<1, 1>>>(aData,textureObj, inputSurfObj);
cudaDeviceSynchronize();
cudaUnbindTexture(&tex);
}
void Aurora::sort(const Aurora::Matrix& aMatrix)
{
RECON_INFO("cuda start");

8
src/Aurora.h
View File

@@ -6,18 +6,10 @@
#include "Matrix.h"
namespace Aurora
{
//__global__ void doubleToComplexKernel(const double* input, cufftDoubleComplex* output, int size);
void doubleToComplex(const double* input, cufftDoubleComplex* output, int size);
//__global__ void maxKernel(const float* aInput, const float* aOutput, int aSize);
void max(const float* aInput, const float* aOutput, int aSize);
Aurora::CudaMatrix valid(const Aurora::CudaMatrix aData, const Aurora::CudaMatrix aValid);
void test(float* aData);
void sort(const Aurora::Matrix& aMatrix);
//Aurora::CudaMatrix getTransmissionDataSubFunction(const Aurora::CudaMatrix& aFxMatrix, const Aurora::CudaMatrix& aFhMatrix);
}

107
src/common/convertfp16tofloat.cpp
View File

@@ -1,107 +0,0 @@
#include "convertfp16tofloat.h"
#include "Function.h"
#include <emmintrin.h>
#include <immintrin.h>
#include <sys/types.h>
namespace {
// const ushort CONVERT_AND_VALUE = 15;
// // andblack
// const __m128i andBlock = _mm_set_epi16(15, 15, 15, 15, 15, 15, 15, 15);
// const __m128i andBlock2 =
// _mm_set_epi16(2047, 2047, 2047, 2047, 2047, 2047, 2047, 2047);
// const __m128i zeroBlock = _mm_set_epi16(0, 0, 0, 0, 0, 0, 0, 0);
// const __m128i oneBlock = _mm_set_epi16(1, 1, 1, 1, 1, 1, 1, 1);
// const __m128i twokBlock =
// _mm_set_epi16(2048, 2048, 2048, 2048, 2048, 2048, 2048, 2048);
// const uint CONVERT_ADD_VALUE = UINT32_MAX - 4095;
// void convert(short * ptr, double* des,bool single = false){
// // 初始化值
// auto value = _mm_set_epi16(ptr[0], ptr[1], ptr[2], ptr[3], single?ptr[0]:ptr[4], single?ptr[0]:ptr[5],
// single?ptr[0]:ptr[6], single?ptr[0]:ptr[7]);
// auto uvalue = _mm_set_epi16(
// (ushort)ptr[0], (ushort)ptr[1], (ushort)ptr[2], (ushort)ptr[3],
// (ushort)(single?ptr[0]:ptr[4]), (ushort)(single?ptr[0]:ptr[5]),
// (ushort)(single?ptr[0]:ptr[6]), (ushort)(single?ptr[0]:ptr[7]));
// // 位移
// auto sign_bit = _mm_srli_epi16(value, 15); // 右移16位取符号位
// auto exponent = _mm_srli_epi16(uvalue, 11);
// // and
// exponent = _mm_and_si128(exponent, andBlock);
// // and ，then convert to int 32 bits
// auto fraction3 = _mm256_cvtepi16_epi32(_mm_and_si128(uvalue, andBlock2));
// auto hidden_bit_mask =
// (_mm_cmp_epi16_mask(sign_bit, oneBlock, _MM_CMPINT_EQ) &
// _mm_cmp_epi16_mask(exponent, zeroBlock, _MM_CMPINT_EQ)) |
// (_mm_cmp_epi16_mask(sign_bit, zeroBlock, _MM_CMPINT_EQ) &
// _mm_cmp_epi16_mask(exponent, zeroBlock, _MM_CMPINT_NE));
// auto hidden_bit16 = _mm_maskz_set1_epi16(hidden_bit_mask, 2048);
// auto hidden_bit32 = _mm256_cvtepi16_epi32(hidden_bit16);
// auto outputBlock = _mm256_add_epi32(fraction3, hidden_bit32);
// auto sign_bit_add_value = _mm256_maskz_set1_epi32(
// _mm_cmp_epi16_mask(sign_bit, oneBlock, _MM_CMPINT_EQ),
// CONVERT_ADD_VALUE);
// outputBlock = _mm256_add_epi32(outputBlock, sign_bit_add_value);
// auto exponent_mask =
// _mm_cmp_epi16_mask(oneBlock, exponent, _MM_CMPINT_LT);
// exponent = _mm_sub_epi16(exponent, oneBlock);
// auto exponent32 = _mm256_cvtepi16_epi32(exponent);
// auto zeroBlock32 = _mm256_cvtepi16_epi32(zeroBlock);
// auto offsetCount =
// _mm256_mask_blend_epi32(exponent_mask, zeroBlock32, exponent32);
// outputBlock = _mm256_sllv_epi32(outputBlock, offsetCount);
// des[3] = _mm256_extract_epi32(outputBlock, 4);
// des[2] = _mm256_extract_epi32(outputBlock, 5);
// des[1] = _mm256_extract_epi32(outputBlock, 6);
// des[0] = _mm256_extract_epi32(outputBlock, 7);
// if(single) return;
// des[7] = _mm256_extract_epi32(outputBlock, 0);
// des[6] = _mm256_extract_epi32(outputBlock, 1);
// des[5] = _mm256_extract_epi32(outputBlock, 2);
// des[4] = _mm256_extract_epi32(outputBlock, 3);
// }
}
Aurora::Matrix Recon::convertfp16tofloat(Aurora::Matrix aMatrix) {
// auto input = aMatrix.getData();
// // uint16变换为float(32位)输出大小翻倍
// auto output = Aurora::malloc(aMatrix.getDataSize() * 4);
// size_t rows = aMatrix.getDataSize() * sizeof(double) / sizeof(short);
// size_t total_count = aMatrix.getDataSize();
// #pragma omp parallel for
// for (size_t i = 0; i < total_count; i += 8) {
// // 循环展开以避免过度的线程调用
// if (i < total_count) {
// auto ptr = (short *)(input + i);
// double *des = output + i * 4;
// ::convert(ptr, des,i+1>total_count);
// }
// if (i+2 < total_count) {
// auto ptr = (short *)(input + i + 2);
// double *des = output + (i+2) * 4;
// ::convert(ptr, des,i+3>total_count);
// }
// if (i+4 < total_count) {
// auto ptr = (short *)(input + i + 4);
// double *des = output + (i+4) * 4;
// ::convert(ptr, des,i+5>total_count);
// }
// if (i+6 < total_count) {
// auto ptr = (short *)(input + i + 6);
// double *des = output + (i+6) * 4;
// ::convert(ptr, des,i+7>total_count);
// }
// }
// return Aurora::Matrix::New(output, aMatrix.getDimSize(0),
// aMatrix.getDimSize(1), aMatrix.getDimSize(2));
}

10
src/common/convertfp16tofloat.h
View File

@@ -1,10 +0,0 @@
#ifndef __CONVERTFP16TOFLOAT_H__
#define __CONVERTFP16TOFLOAT_H__
#include "Matrix.h"
namespace Recon {
Aurora::Matrix convertfp16tofloat(Aurora::Matrix aMatrix);
}
#endif // __CONVERTFP16TOFLOAT_H__

98
src/common/dataBlockCreation/getAScanBlockPreprocessed.cpp
View File

@@ -1,6 +1,5 @@
#include "getAScanBlockPreprocessed.h"
#include "CudaMatrix.h"
#include "Matrix.h"
#include "blockingGeometryInfo.h"
#include "removeDataFromArrays.h"
@@ -11,36 +10,15 @@
#include "src/transmissionReconstruction/dataFilter/dataFilter.h"
#include "src/reflectionReconstruction/dataFilter.h"
#include "Aurora.h"
using namespace Aurora;
using namespace Recon;
#include <sys/time.h>
#include <iostream>
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;
}
AscanBlockPreprocessed Recon::getAscanBlockPreprocessed(Parser* aParser, const Aurora::Matrix& aMp, const Aurora::Matrix& aSl, const Aurora::Matrix& aSn,
const Aurora::Matrix& aRl, const Aurora::Matrix& aRn, GeometryInfo& aGeom, const MeasurementInfo& aMeasInfo,
bool aApplyFilter, bool aTransReco)
{
//std::cout<<"strart"<<std::endl;
//printTime();
//550ms
AscanBlockPreprocessed result;
AscanBlock ascanBlock = getAscanBlock(aParser, aMp, aSl, aSn, aRl, aRn);
//printTime();
//10ms
result.gainBlock = ascanBlock.gainBlock;
result.mpBlock = ascanBlock.mpBlock;
result.rlBlock = ascanBlock.rlBlock;
@@ -48,8 +26,6 @@ AscanBlockPreprocessed Recon::getAscanBlockPreprocessed(Parser* aParser, const A
result.slBlock = ascanBlock.slBlock;
result.snBlock = ascanBlock.snBlock;
GeometryBlock geometryBlock = blockingGeometryInfos(aGeom, ascanBlock.rnBlock, ascanBlock.rlBlock, ascanBlock.snBlock, ascanBlock.slBlock, ascanBlock.mpBlock);
//printTime();
//3ms
result.receiverPositionBlock = geometryBlock.receiverPositionBlock;
result.senderPositionBlock = geometryBlock.senderPositionBlock;
if(aApplyFilter)
@@ -64,8 +40,7 @@ AscanBlockPreprocessed Recon::getAscanBlockPreprocessed(Parser* aParser, const A
{
usedData = filterReflectionData(geometryBlock.receiverPositionBlock, geometryBlock.senderPositionBlock, geometryBlock.senderNormalBlock, reflectParams::constrictReflectionAngles);
}
//printTime();
//150ms
ascanBlock.ascanBlock = removeDataFromArrays(ascanBlock.ascanBlock, usedData);
result.mpBlock = removeDataFromArrays(ascanBlock.mpBlock, usedData);
result.slBlock = removeDataFromArrays(ascanBlock.slBlock, usedData);
@@ -76,8 +51,7 @@ AscanBlockPreprocessed Recon::getAscanBlockPreprocessed(Parser* aParser, const A
result.senderPositionBlock = removeDataFromArrays(geometryBlock.senderPositionBlock, usedData);
result.receiverPositionBlock = removeDataFromArrays(geometryBlock.receiverPositionBlock, usedData);
result.gainBlock = removeDataFromArrays(ascanBlock.gainBlock, usedData);
//printTime();
//120ms
}
if (ascanBlock.ascanBlock.getDataSize() > 0)
@@ -88,72 +62,6 @@ AscanBlockPreprocessed Recon::getAscanBlockPreprocessed(Parser* aParser, const A
{
result.ascanBlockPreprocessed = ascanBlock.ascanBlock;
}
//printTime();
return result;
}
AscanBlockPreprocessedCuda Recon::getAscanBlockPreprocessedCuda(Parser* aParser, const Aurora::Matrix& aMp, const Aurora::Matrix& aSl, const Aurora::Matrix& aSn,
const Aurora::Matrix& aRl, const Aurora::Matrix& aRn, GeometryInfo& aGeom, const MeasurementInfo& aMeasInfo,
bool aApplyFilter, bool aTransReco)
{
//std::cout<<"strart"<<std::endl;
//printTime();
//550ms
AscanBlockPreprocessedCuda result;
AscanBlock ascanBlock = getAscanBlock(aParser, aMp, aSl, aSn, aRl, aRn);
//printTime();
//300ms
result.ascanBlockPreprocessed = ascanBlock.ascanBlock.toDeviceMatrix();
result.gainBlock = ascanBlock.gainBlock.toDeviceMatrix();
result.mpBlock = ascanBlock.mpBlock;
result.rlBlock = ascanBlock.rlBlock;
result.rnBlock = ascanBlock.rnBlock;
result.slBlock = ascanBlock.slBlock;
result.snBlock = ascanBlock.snBlock;
GeometryBlock geometryBlock = blockingGeometryInfos(aGeom, ascanBlock.rnBlock, ascanBlock.rlBlock, ascanBlock.snBlock, ascanBlock.slBlock, ascanBlock.mpBlock);
//printTime();
//3ms
result.receiverPositionBlock = geometryBlock.receiverPositionBlock;
result.senderPositionBlock = geometryBlock.senderPositionBlock;
if(aApplyFilter)
{
Matrix usedData;
if(aTransReco)
{
usedData = filterTransmissionData(ascanBlock.slBlock, ascanBlock.snBlock, ascanBlock.rlBlock, ascanBlock.rnBlock,
aGeom.sensData, geometryBlock.senderNormalBlock, geometryBlock.receiverNormalBlock);
}
else
{
usedData = filterReflectionData(geometryBlock.receiverPositionBlock, geometryBlock.senderPositionBlock, geometryBlock.senderNormalBlock, reflectParams::constrictReflectionAngles);
}
//printTime();
//40ms
CudaMatrix usedDataDevice = usedData.toDeviceMatrix();
result.ascanBlockPreprocessed = valid(result.ascanBlockPreprocessed, usedDataDevice);
result.mpBlock = removeDataFromArrays(ascanBlock.mpBlock, usedData);
result.slBlock = removeDataFromArrays(ascanBlock.slBlock, usedData);
result.snBlock = removeDataFromArrays(ascanBlock.snBlock, usedData);
result.rlBlock = removeDataFromArrays(ascanBlock.rlBlock, usedData);
result.rnBlock = removeDataFromArrays(ascanBlock.rnBlock, usedData);
result.senderPositionBlock = removeDataFromArrays(geometryBlock.senderPositionBlock, usedData);
result.receiverPositionBlock = removeDataFromArrays(geometryBlock.receiverPositionBlock, usedData);
result.gainBlock = valid(result.gainBlock, usedDataDevice);
//printTime();
//10ms
}
if (ascanBlock.ascanBlock.getDataSize() > 0)
{
result.ascanBlockPreprocessed = preprocessAscanBlockCuda(result.ascanBlockPreprocessed, aMeasInfo);
}
// else
// {
// result.ascanBlockPreprocessed = ascanBlock.ascanBlock;
// }
//printTime();
//std::cout<<"end"<<std::endl;
return result;
}

18
src/common/dataBlockCreation/getAScanBlockPreprocessed.h
View File

@@ -2,7 +2,6 @@
#define GETASCANBLOCK_PREPROCESSED_H
#include "Matrix.h"
#include "CudaMatrix.h"
#include "src/common/getGeometryInfo.h"
#include "src/common/getMeasurementMetaData.h"
@@ -26,23 +25,6 @@ namespace Recon
AscanBlockPreprocessed getAscanBlockPreprocessed(Parser* aParser, const Aurora::Matrix& aMp, const Aurora::Matrix& aSl, const Aurora::Matrix& aSn,
const Aurora::Matrix& aRl, const Aurora::Matrix& aRn, GeometryInfo& aGeom, const MeasurementInfo& aMeasInfo,
bool aApplyFilter, bool aTransReco);
struct AscanBlockPreprocessedCuda
{
Aurora::CudaMatrix ascanBlockPreprocessed;
Aurora::Matrix mpBlock;
Aurora::Matrix slBlock;
Aurora::Matrix snBlock;
Aurora::Matrix rlBlock;
Aurora::Matrix rnBlock;
Aurora::Matrix senderPositionBlock;
Aurora::Matrix receiverPositionBlock;
Aurora::CudaMatrix gainBlock;
};
AscanBlockPreprocessedCuda getAscanBlockPreprocessedCuda(Parser* aParser, const Aurora::Matrix& aMp, const Aurora::Matrix& aSl, const Aurora::Matrix& aSn,
const Aurora::Matrix& aRl, const Aurora::Matrix& aRn, GeometryInfo& aGeom, const MeasurementInfo& aMeasInfo,
bool aApplyFilter, bool aTransReco);
}

4
src/common/fileHelper.h
View File

@@ -6,8 +6,8 @@
namespace Recon
{
const std::string DEFAULT_CONFIG_PATH = "/home/UR/ConfigFiles/";
const std::string DEFAULT_OUTPUT_PATH = "/home/UR/ReconResult/";
const std::string DEFAULT_OUTPUT_FILENAME = "sun.mat";
const std::string DEFAULT_OUTPUT_PATH = "/home/UR/ReconResult/USCT_Result.mat";
const std::string DEFAULT_OUTPUT_FILENAME = "USCT_Result.mat";
std::string getPath(const std::string &aFullPath);
bool endsWithMat(const std::string &aStr);

16
src/log/notify.cpp
View File

@@ -1,16 +0,0 @@
#include "notify.h"
#include "log/log.h"
namespace Recon {
std::string ReconID;
bool notifyStart(const std::string& aReconID ){
return true;
}
bool notifyFinish(){
return true;
}
bool notifyProgress( int percent){
return true;
}
}

9
src/log/notify.h
View File

@@ -1,9 +0,0 @@
#ifndef __NOTIFY_H__
#define __NOTIFY_H__
#include <string>
namespace Recon {
bool notifyStart(const std::string& aReconID );
bool notifyFinish();
bool notifyProgress(int percent);
}
#endif // __NOTIFY_H__

13
src/main.cxx
View File

@@ -1,19 +1,12 @@
#include "Parser.h"
#include "Parser.h"
#include "config/config.h"
#include "log/notify.h"
#include "log/notify.h"
#include <cstdio>
#include <vector>
#include "MatlabReader.h"
#define EIGEN_USE_MKL_ALL
#include "src/common/fileHelper.h"
#include "startReconstructions.h"
#include "transmissionReconstruction/detection/detection.h"
#include "transmissionReconstruction/detection/detection.cuh"
#include "startReconstructions.h"
#include "log/log.h"
/* 0 is data path.
1 is dataRef path.
@@ -25,9 +18,9 @@ int main(int argc, char *argv[])
{
int argNum = 5;
std::vector<std::string> args(argNum);
args[0] = "/home/AScans_Data/CAS0.1Phantom/20230823T165617/";
args[1] = "/home/AScans_Data/CAS0.1Phantom/20230823T171851/";
args[2] = "/home/krad/Storage/DICOM/00e04b741e9f_20240619T145752/";
args[0] = "/home/sun/20230418T145123/";
args[1] = "/home/sun/20230418T141000/";
args[2] = "/home/sun/AscanData/";
args[3] = Recon::DEFAULT_CONFIG_PATH;
argc = argc <= argNum? argc-1 : argNum;
for (int i = 0; i < argc; i++)

2
src/reflectionReconstruction/startReflectionReconstruction.cpp
View File

@@ -9,7 +9,6 @@
#include "common/precalculateChannelList.h"
#include "common/dataBlockCreation/getAScanBlockPreprocessed.h"
#include "common/dataBlockCreation/removeDataFromArrays.h"
#include "log/notify.h"
#include "reflectionReconstruction/preprocessData/determineOptimalPulse.h"
#include "reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.h"
#include "src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.h"
@@ -139,7 +138,6 @@ Aurora::Matrix Recon::startReflectionReconstruction( Parser* aParser, int aSAFT_
std::cout<<Env[0]<<"-" << Env[1] <<"-" << Env[2] <<"-" << Env[3]<<std::endl;
RECON_INFO("Reflection Reconstructon: " + std::to_string(j));
}
Recon::notifyProgress(25+73*((j*i)/(aMotorPos.getDataSize() * aSlList.getDataSize())));
}
}

20
src/startReconstructions.cpp
View File

@@ -8,7 +8,6 @@
#include "common/getMeasurementMetaData.h"
#include "common/getGeometryInfo.h"
#include "common/estimatePulseLength.h"
#include "log/notify.h"
#include "transmissionReconstruction/startTransmissionReconstruction.h"
#include <string>
@@ -33,7 +32,7 @@ using namespace Aurora;
int Recon::startReconstructions( const std::string& aDataPath, const std::string& aDataRefPath, const std::string& aOutputPath)
{
// MatlabWriter writer(aOutputPath);
MatlabWriter writer(aOutputPath + "URResult.mat");
Parser dataParser(aDataPath);
Parser refParser(aDataRefPath);
Recon::DICOMExporter exporter(dataParser.getPatientData(), dataParser.getMetaData());
@@ -96,7 +95,6 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
TempInfo tempRef;
CEInfo ceRef;
Matrix transformationMatricesRef;
//Recon::notifyProgress(1);
if(transParams::runTransmissionReco)
{
expInfoRef = loadMeasurementInfos(&refParser);
@@ -129,7 +127,6 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
transformationMatricesRef = Matrix();
motorPosAvailableRef = Matrix();
}
//Recon::notifyProgress(2);
if(!ce.ce.isNull() && !ceRef.ce.isNull())
{
Matrix isEqual = (ce.ce == ceRef.ce);
@@ -156,12 +153,10 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
preComputes.matchedFilter = createMatchedFilter(ce.ceRef, ce.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
}
}
//Recon::notifyProgress(3);
if(expInfo.sampleRate != reflectParams::aScanReconstructionFrequency)
{
reflectParams::expectedAScanDataLength = ceil(expInfo.numberSamples * ((float)reflectParams::aScanReconstructionFrequency / expInfo.sampleRate));
}
//Recon::notifyProgress(4);
TransmissionReconstructionResult transmissionResult;
bool sosAvailable = false;
bool attAvailable = false;
@@ -220,14 +215,13 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
GeometryInfo geomRef = getGeometryInfo(motorPosAvailableRef, transformationMatricesRef, rlList, rnList, slList, snList);
RECON_INFO("Start transmissionRecostruction.");
//Recon::notifyProgress(5);
transmissionResult = startTransmissionReconstruction(mp_inter, mpRef_inter, slList_inter, snList_inter, rlList_inter, rnList_inter, temp, tempRef, geom, geomRef, expInfo, expInfoRef, preComputes, &dataParser, &refParser);
attAvailable = true;
sosAvailable = true;
exporter.exportDICOM(transmissionResult.recoSOS, DICOMExporter::SOS);
exporter.exportDICOM(transmissionResult.recoATT, DICOMExporter::ATT);
// writer.setMatrix(transmissionResult.recoSOS, "sos");
// writer.setMatrix(transmissionResult.recoATT, "att");
writer.setMatrix(transmissionResult.recoSOS, "sos");
writer.setMatrix(transmissionResult.recoATT, "att");
}
if(reflectParams::runReflectionReco)
@@ -235,11 +229,9 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
Matrix recoSOS = transmissionResult.recoSOS;
Matrix recoATT = transmissionResult.recoATT;
precalcImageParameters(geom);
//Recon::notifyProgress(21);
//检测可使用内存是否足够,输出警报用,todo
//checkEnvAndMemory(reflectParams.imageInfos.IMAGE_XYZ);
auto preProcessData = preprocessTransmissionReconstructionForReflection(recoSOS, recoATT, transmissionResult.ddmis, geom, temp);
//Recon::notifyProgress(22);
Matrix mp_inter = intersect(motorPosAvailable, reflectParams::motorPos);
Matrix slList_inter = intersect(slList, reflectParams::senderTasList);
Matrix snList_inter = intersect(snList, reflectParams::senderElementList);
@@ -252,14 +244,12 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
preComputes.offset = estimateOffset(expInfo, ce, preComputes.matchedFilter);
reflectParams::gpuSelectionList = reconParams::gpuSelectionList;
//Recon::notifyProgress(25);
RECON_INFO("Start reflectionRecostruction.");
Matrix env = startReflectionReconstruction(&dataParser, preProcessData.saftMode, mp_inter, slList_inter, snList_inter, rlList_inter, rnList_inter, geom, preProcessData.transRecos, expInfo, preComputes);
// writer.setMatrix(env, "reflect");
writer.setMatrix(env, "reflect");
exporter.exportDICOM(env, Recon::DICOMExporter::REFL);
//Recon::notifyProgress(99);
}
// writer.write();
writer.write();
return 0;
}

3
src/transmissionReconstruction/dataPreperation.h
View File

@@ -3,8 +3,7 @@
#include "Matrix.h"
namespace Recon {
Aurora::Matrix distanceBetweenTwoPoints(const Aurora::Matrix& aMPtsA,
const Aurora::Matrix& aMPtsB);
Aurora::Matrix distanceBetweenTwoPoints(const Aurora::Matrix& aMPtsA, const Aurora::Matrix& aMPtsB);
Aurora::Matrix calculateWaterTemperature(Aurora::Matrix aMWaterTempS,
Aurora::Matrix aMWaterTempR,

42
src/transmissionReconstruction/detection/detection.cpp
View File

@@ -13,7 +13,6 @@
#include "common/getMeasurementMetaData.h"
#include "config/config.h"
#include "calculateBankDetectAndHilbertTransformation.hpp"
#include "transmissionReconstruction/detection/detection.cuh"
using namespace Aurora;
namespace Recon {
@@ -415,39 +414,34 @@ namespace Recon {
return result;
}
DetectResult transmissionDetection(const Aurora::CudaMatrix &AscanBlock,
const Aurora::CudaMatrix &AscanRefBlock,
const Aurora::CudaMatrix &distBlock,
const Aurora::CudaMatrix &distRefBlock,
DetectResult transmissionDetection(const Aurora::Matrix &AscanBlock,
const Aurora::Matrix &AscanRefBlock,
const Aurora::Matrix &distBlock,
const Aurora::Matrix &distRefBlock,
const Aurora::Matrix &sosWaterBlock,
const Aurora::Matrix &sosWaterRefBlock,
float expectedSOSWater) {
auto _sosWaterBlock = temperatureToSoundSpeed(sosWaterBlock, "marczak").toDeviceMatrix();
auto _sosWaterRefBlock = temperatureToSoundSpeed(sosWaterRefBlock, "marczak").toDeviceMatrix();
auto _sosWaterBlock = temperatureToSoundSpeed(sosWaterBlock, "marczak");
auto _sosWaterRefBlock = temperatureToSoundSpeed(sosWaterRefBlock, "marczak");
switch (Recon::transParams::version) {
// case 1: {
// return detectTofAndAttMex(
// AscanBlock, AscanRefBlock, distBlock, distRefBlock,
// _sosWaterBlock, _sosWaterRefBlock, Recon::transParams::resampleFactor, Recon::transParams::nThreads,
// expectedSOSWater, Recon::transParams::useTimeWindowing,
// Recon::transParams::aScanReconstructionFrequency, Recon::transParams::detectionWindowATT,
// Recon::transParams::offsetElectronic, Recon::transParams::detectionWindowSOS, Recon::transParams::minSpeedOfSound,
// Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow);
// }
// case 2:
default:
auto r = detectTofAndAtt(
case 1: {
return detectTofAndAttMex(
AscanBlock, AscanRefBlock, distBlock, distRefBlock,
_sosWaterBlock, _sosWaterRefBlock, Recon::transParams::resampleFactor, Recon::transParams::nThreads,
expectedSOSWater, Recon::transParams::useTimeWindowing,
Recon::transParams::aScanReconstructionFrequency, Recon::transParams::detectionWindowATT,
Recon::transParams::offsetElectronic, Recon::transParams::detectionWindowSOS, Recon::transParams::minSpeedOfSound,
Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow);
}
case 2:
default:
return detectTofAndAtt(
AscanBlock, AscanRefBlock, distBlock, distRefBlock,
_sosWaterBlock, _sosWaterRefBlock, Recon::transParams::resampleFactor, Recon::transParams::nThreads,
expectedSOSWater, Recon::transParams::useTimeWindowing,
Recon::transParams::aScanReconstructionFrequency, Recon::transParams::detectionWindowATT,
Recon::transParams::offsetElectronic, Recon::transParams::detectionWindowSOS, Recon::transParams::minSpeedOfSound,
Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow);
DetectResult ret;
ret.att = r.att.toHostMatrix();
ret.sosValue = r.sosValue.toHostMatrix();
ret.tof = r.tof.toHostMatrix();
return ret;
}
}
}

4
src/transmissionReconstruction/detection/detection.h
View File

@@ -75,8 +75,8 @@ detectTofAndAttMex(
DetectResult
transmissionDetection(
const Aurora::CudaMatrix &AscanBlock, const Aurora::CudaMatrix &AscanRefBlock,
const Aurora::CudaMatrix &distBlock, const Aurora::CudaMatrix &distRefBlock,
const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
const Aurora::Matrix &distBlock, const Aurora::Matrix &distRefBlock,
const Aurora::Matrix &sosWaterBlock, const Aurora::Matrix &sosWaterRefBlock, float expectedSOSWater);
} // namespace Recon

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

@@ -1,12 +1,8 @@
#include "getTransmissionData.h"
#include "getTransmissionData.cuh"
#include "AuroraDefs.h"
#include "CudaMatrix.h"
#include "Function.h"
#include "Function1D.h"
#include "Function2D.h"
#include "common/dataBlockCreation/removeDataFromArrays.h"
#include "log/notify.h"
#include "src//config/config.h"
#include "src/common/getGeometryInfo.h"
#include "src/common/temperatureCalculation/extractTasTemperature.h"
@@ -31,11 +27,6 @@
#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;
@@ -50,8 +41,8 @@ namespace
Matrix waterTempBlock;
MetaInfos metaInfos;
CudaMatrix ascanBlock;
CudaMatrix ascanBlockRef;
Matrix ascanBlock;
Matrix ascanBlockRef;
Matrix dists;
Matrix distRefBlock;
Matrix waterTempRefBlock;
@@ -66,32 +57,13 @@ namespace
int BUFFER_COUNT = 0;
int BUFFER_SIZE = 4;//<=8
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 prepareAScansForTransmissionDetection(const Matrix& aAscanBlock, const Matrix& aGainBlock)
{
CudaMatrix result = aAscanBlock / repmat(aGainBlock, aAscanBlock.getDimSize(0), 1);
Matrix 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,
@@ -99,24 +71,12 @@ void printTime()
{
BlockOfTransmissionData result;
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();
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();
if(aExpInfo.Hardware == "USCT3dv3")
{
Matrix channelList = precalculateChannelList(aRlList, aRnList, aExpInfo, aPreComputes);
@@ -132,48 +92,71 @@ void printTime()
channelListBlockData[i] = channelList[ind[i] - 1];
}
Matrix channelListBlock = Matrix::New(channelListBlockData, 1, channelListBlockSize);
//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);
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);
ascanBlock = Aurora::real(ifft(complex));
//printTime();
//20s
fx = fft(ascanBlockRef);
//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);
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);
ascanBlockRef = Aurora::real(ifft(complex));
//printTime();
//20ms
}
else
{
@@ -183,27 +166,24 @@ void printTime()
if(transParams::applyCalib)
{
metaInfos.snrValues = calculateSnr(ascanBlock, aSnrRmsNoise[0]).toHostMatrix();
metaInfos.snrValuesRef = calculateSnr(ascanBlockRef, aSnrRmsNoiseRef[0]).toHostMatrix();
metaInfos.snrValues = calculateSnr(ascanBlock, aSnrRmsNoise[0]);
metaInfos.snrValuesRef = calculateSnr(ascanBlockRef, aSnrRmsNoiseRef[0]);
}
// 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);
// 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;
// }
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;
@@ -218,7 +198,7 @@ void printTime()
// DetectResult detect = transmissionDetection(ascanBlock, ascanBlockRef, dists, distRefBlock, waterTempBlock, waterTempRefBlock, aExpectedSOSWater[0]);
// result.attData = detect.att;
// result.tofData = detect.tof;
//printTime();
return result;
}
@@ -227,9 +207,8 @@ void printTime()
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, unsigned int aGPUId)
const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
{
cudaSetDevice(aGPUId);
auto buffer = getBlockOfTransmissionData(aMp, aMpRef, aSl, aSn, aRlList, aRnList, aTasTemps,
aExpectedSOSWater, aGeom, aGeomRef, aSnrRmsNoise, aSnrRmsNoiseRef,
aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef);
@@ -265,7 +244,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, index % BUFFER_SIZE);
speedUpThread[index] = std::thread(getBlockOfTransmissionDataInThread,index,mp,mpRef,sl,sn,aRlList,aRnList,aTasTemps,aExpectedSOSWater,aGeom,aGeomRef,aSnrRmsNoise,aSnrRmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef);
}
}
}
@@ -354,9 +333,8 @@ 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.toDeviceMatrix(), blockData.distRefBlock.toDeviceMatrix(),
blockData.dists, blockData.distRefBlock,
blockData.waterTempBlock, blockData.waterTempRefBlock,
aTemp.expectedSOSWater[0]);
blockData.attData = detect.att;
@@ -376,14 +354,6 @@ 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);
// }
numData += numUsedData;
std::unique_lock<std::mutex> lockBufferCount(CREATE_BUFFER_MUTEX);
BLOCK_OF_TRANSIMISSIONDARA_BUFFER.erase(std::to_string(index));
@@ -392,7 +362,6 @@ 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))));
}
}
speedUpThread.join();
@@ -420,14 +389,6 @@ 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);
// }
Matrix valid;
if(transParams::applyCalib)
@@ -453,15 +414,6 @@ 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;
// }
tofDataTotal = removeDataFromArrays(tofDataTotal, valid);
attDataTotal = removeDataFromArrays(attDataTotal, valid);
senderList = removeDataFromArrays(senderList, valid);

27
src/transmissionReconstruction/reconstruction/reconstruction.cpp
View File

@@ -7,21 +7,14 @@
#include "config/config.h"
#include "CudaEnvInit.h"
#include "log/notify.h"
#include "transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h"
#include "transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cuh"
#include "src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.h"
#include <algorithm>
#include <cmath>
#include <iostream>
#include <vector>
#include <thread>
#include <future>
using namespace Aurora;
using solveParameterIteratorFunctionType = std::vector<std::vector<Aurora::Matrix>> (*)(Aurora::Sparse M, Aurora::Matrix &b,
const Aurora::Matrix &dims, bool oneIter, bool nonNeg, int aDevice);
using slownessToSOSFunctionType = Matrix (*)(Aurora::Matrix & aVF1, float aSOS_IN_WATER);
namespace Recon {
Aurora::Matrix calculateMinimalMaximalTransducerPositions(
const Aurora::Matrix &aMSenderList, const Aurora::Matrix &aMReceiverList) {
@@ -195,9 +188,7 @@ namespace Recon {
BuildMatrixResult buildMatrixR;
for(int iter=1; iter<=numIter; ++iter)
{
auto resDevice = res.toDeviceMatrix();
//1200ms
buildMatrixR = buildMatrix(senderList.toDeviceMatrix(), receiverList.toDeviceMatrix(), resDevice, dims.toDeviceMatrix(), bentRecon && (iter!=1), potentialMap.toDeviceMatrix());
buildMatrixR = buildMatrix(senderList, receiverList, res, dims, bentRecon && (iter!=1), potentialMap);
if(!data.isNull() && bentRecon && iter != numIter)
{
//与默认配置bentRecon不符，暂不实现todo
@@ -230,23 +221,25 @@ namespace Recon {
{
allHitMaps.push_back(buildMatrixR.hitmap);
}
#pragma omp parallel for num_threads(2)
for (int i =0; i<2; i++){
if (i ==0){
if(!data.isNull())
{
//1500ms
Matrix sosValue = solveParameterIterator(buildMatrixR.M, b, dims, false, transParams::nonNeg)[0][0];
//1ms
result.outSOS = slownessToSOS(sosValue, SOS_IN_WATER) ;
}
}
else{
if(!dataAtt.isNull())
{
//1500ms
Matrix attValue = solveParameterIterator(buildMatrixR.M, bAtt, dims, false, transParams::nonNeg)[0][0];
//1ms
Matrix attValue = solveParameterIterator(buildMatrixR.M, bAtt, dims, false, transParams::nonNeg,1)[0][0];
result.outATT = attValue/100 ;
}
}
}
}
return result;
}

6
src/transmissionReconstruction/startTransmissionReconstruction.cpp
View File

@@ -1,10 +1,8 @@
#include "startTransmissionReconstruction.h"
#include "./detection/getTransmissionData.h"
#include "Matrix.h"
#include "CudaMatrix.h"
#include "log/log.h"
#include "common/dataBlockCreation/removeDataFromArrays.h"
#include "log/notify.h"
#include "src/transmissionReconstruction/dataFilter/dataFilter.h"
#include "src/transmissionReconstruction/dataPreperation.h"
#include "src/common/getMeasurementMetaData.h"
@@ -30,22 +28,18 @@ TransmissionReconstructionResult Recon::startTransmissionReconstruction(const Au
aGeom, aGeomRef, aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef);
Matrix dists = Recon::distanceBetweenTwoPoints(transmissionData.senderList, transmissionData.receiverList);
Matrix sosRef = Recon::temperatureToSoundSpeed(transmissionData.waterTempList, "marczak");
//Recon::notifyProgress(17);
Matrix valid = Recon::checkTofDetections(transmissionData.tofDataTotal, dists, sosRef,
Recon::transParams::minSpeedOfSound,Recon::transParams::maxSpeedOfSound).valid;
//Recon::notifyProgress(18);
if(transParams::qualityCheck)
{
qualityReview(sum(valid,Aurora::All)[0], transmissionData.dataInfo.numPossibleScans);
}
//Recon::notifyProgress(19);
DiscretizePositionValues positionValues = Recon::discretizePositions(transmissionData.senderList, transmissionData.receiverList, Recon::transParams::numPixelXY);
Matrix tofData = removeDataFromArrays(transmissionData.tofDataTotal, valid);
Matrix attData = removeDataFromArrays(transmissionData.attDataTotal, valid);
Matrix senderList = removeDataFromArrays(positionValues.senderCoordList, valid);
Matrix reveiverList = removeDataFromArrays(positionValues.receiverCoordList, valid);
//Recon::notifyProgress(20);
RECON_INFO("Start reconstructArt.");
auto transmissionReon = reconstructArt(tofData, attData, positionValues.dims, senderList, reveiverList, positionValues.res, aTemp.expectedSOSWater[0]);

16
test/Common_Test.cpp
View File

@@ -4,7 +4,6 @@
#include "Matrix.h"
#include "common/ceMatchedFilterHandling.h"
#include "common/common.h"
#include "common/convertfp16tofloat.h"
#include "common/getGeometryInfo.h"
#include "common/dataBlockCreation/getAscanBlock.h"
#include "common/dataBlockCreation/blockingGeometryInfo.h"
@@ -52,21 +51,6 @@ TEST_F(Common_Test, adaptFrequency) {
}
TEST_F(Common_Test, convertfp16tofloat) {
MatlabReader m("/home/krad/TestData/convertReal.mat");
size_t count = 0;
auto input = m.readint16("input",count);
auto ma = Aurora::Matrix::copyFromRawData((float*)input.get(),count/4);
auto resultM = Recon::convertfp16tofloat(ma);
auto result = resultM.getData();
auto output = m.read("output");
for (size_t i = 0; i<count; i++) {
EXPECT_EQ(result[i], output.getData()[i])<<"index:"<<i<<",input:"<< ((short*)ma.getData())[i]<<",input2:"<<input.get()[i];
}
}
TEST_F(Common_Test, convertToLinearIndices) {
//2D
{