diff --git a/CMakeLists.txt b/CMakeLists.txt index 30033d5..e38ab0e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,7 +6,7 @@ set(CMAKE_CXX_STANDARD 14) set(CMAKE_INCLUDE_CURRENT_DIR ON) set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -g") - +set(MKL_LINK static) find_package(Aurora REQUIRED) find_package(Parser REQUIRED) file(GLOB_RECURSE cpp_files ./src/*.cpp) @@ -22,4 +22,5 @@ target_include_directories(UR PUBLIC ${Parser_INCLUDE_DIRS}) target_link_libraries(UR PUBLIC ${Aurora_Libraries}) target_link_libraries(UR PUBLIC matio) +target_link_libraries(UR PUBLIC pthread) target_link_libraries(UR PUBLIC Parser) \ No newline at end of file diff --git a/src/common/DICOMExporter.cpp b/src/common/DICOMExporter.cpp deleted file mode 100644 index f608dc7..0000000 --- a/src/common/DICOMExporter.cpp +++ /dev/null @@ -1,260 +0,0 @@ -#include "DICOMExporter.h" - -#include "fileHelper.h" -#include "config/config.h" -#include "log/log.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include "dcmtk/dcmdata/dctk.h" - -using namespace std; -namespace { - float getSpacing(float * startPoint, float* endPoint, int* imageXYZ, int index){ - return (endPoint[index]-startPoint[index])/(float)(imageXYZ[index])*1000.0; -} - - void initialDicomFile(DcmDataset* dataset,DcmMetaInfo* metaInfo) - { - metaInfo->putAndInsertString(DCM_MediaStorageSOPClassUID, - UID_CTImageStorage); - dataset->putAndInsertString(DCM_SOPClassUID, UID_CTImageStorage); - - dataset->putAndInsertString(DCM_Modality, "CT"); - dataset->putAndInsertString(DCM_TransferSyntaxUID, - UID_LittleEndianImplicitTransferSyntax); - dataset->putAndInsertString(DCM_ImageType, "ORIGINAL\\PRIMARY"); - - dataset->putAndInsertString(DCM_Manufacturer, "EQ9"); - dataset->putAndInsertString(DCM_PatientOrientation, "L\\H"); - dataset->putAndInsertString(DCM_PatientPosition, "HFP"); - - dataset->putAndInsertString(DCM_ImageOrientationPatient, "1\\0\\0\\0\\0\\1"); - - dataset->putAndInsertString(DCM_StudyDescription, "USCT BREAST"); - dataset->putAndInsertString(DCM_ManufacturerModelVersion, "USCT3DV3"); - - dataset->putAndInsertString(DCM_SamplesPerPixel, "1"); - dataset->putAndInsertString(DCM_PixelRepresentation, "0"); - dataset->putAndInsertString(DCM_BitsAllocated, "16"); - dataset->putAndInsertString(DCM_BitsStored, "16"); - dataset->putAndInsertString(DCM_HighBit, "15"); - dataset->putAndInsertString(DCM_PhotometricInterpretation, "MONOCHROME2"); - - dataset->putAndInsertString(DCM_KVP, "0"); - - } - void initialDate(DcmDataset* dataset,const std::string& mStudyDate ) - { - dataset->putAndInsertString(DCM_StudyDate, mStudyDate.data()); - OFDateTime currentDateTime; - currentDateTime.setCurrentDateTime(); - OFString date; - currentDateTime.getDate().getISOFormattedDate(date,false); - dataset->putAndInsertString(DCM_SeriesDate, date.data()); - dataset->putAndInsertString(DCM_AcquisitionDate, date.data()); - dataset->putAndInsertString(DCM_ContentDate, date.data()); - } - - void initialSeriesTime(DcmDataset* dataset){ - // 获取当前日期和时间(UTC时间) - OFDateTime currentDateTime; - currentDateTime.setCurrentDateTime(); - OFString time; - currentDateTime.getTime().getISOFormattedTime(time,true,false,false,false); - dataset->putAndInsertString(DCM_SeriesTime, time.data()); - dataset->putAndInsertString(DCM_AcquisitionTime, time.data()); - dataset->putAndInsertString(DCM_ContentTime, time.data()); - } - const char* SeriesDescriptions[3]={"REFL","SOS","ATT"}; - const char* getSeriesDescription(int type){ - return SeriesDescriptions[type]; - } -} - - -namespace Recon -{ - DICOMExporter::DICOMExporter(const PatientData& aPatientData) - { - mPatientData = aPatientData; - //prepare Study - dcmGenerateUniqueIdentifier(mStudyInstanceUID,SITE_STUDY_UID_ROOT); - - OFDateTime currentDateTime; - currentDateTime.setCurrentDateTime(); - OFString date; - currentDateTime.getDate().getISOFormattedDate(date,false); - OFString time; - currentDateTime.getTime().getISOFormattedTime(time,true,false,false,false); - mStudyDate = date.data(); - mStudyTime = time.data(); - } - - DICOMExporter::DICOMExporter() - { - //prepare Study - dcmGenerateUniqueIdentifier(mStudyInstanceUID,SITE_STUDY_UID_ROOT); - - OFDateTime currentDateTime; - currentDateTime.setCurrentDateTime(); - OFString date; - currentDateTime.getDate().getISOFormattedDate(date,false); - OFString time; - currentDateTime.getTime().getISOFormattedTime(time,true,false,false,false); - mStudyDate = date.data(); - mStudyTime = time.data(); - } - - DICOMExporter::~DICOMExporter() - { - } - void DICOMExporter::setExportBasePath(const std::string & path) - { - mBasePath = path; - } - - void DICOMExporter::exportDICOM(Aurora::Matrix aMatrix, ImageType type){ - int XYZ [3] = {aMatrix.getDimSize(0),aMatrix.getDimSize(1), aMatrix.getDimSize(2)}; - string path = mBasePath+"/"+::getSeriesDescription(type)+"/"; - if (!exists(path)){ - if(!mkdir(path)){ - RECON_ERROR("Fail to create dir with path {0}\r\n",path); - } - } - exportDCM(path, reflectParams::imageStartpoint.getData(), - reflectParams::imageEndpoint.getData(), XYZ, aMatrix.getData(), type); - RECON_INFO("Save DICOM to Path:{0}, type:{1}==================>",path,(int)type); - } - - void DICOMExporter::exportDCM(string path, float * startPoint, float* endPoint, int* imageXYZ,float* data, int type) - { - long Rows = imageXYZ[1], Cols = imageXYZ[0] ,Slices = imageXYZ[2]; - DcmFileFormat dcmFile; - DcmDataset* dataset = dcmFile.getDataset(); - DcmMetaInfo* metaInfo = dcmFile.getMetaInfo(); - - initialDicomFile(dataset,metaInfo); - - dataset->putAndInsertString(DCM_StudyInstanceUID, mStudyInstanceUID); - - dcmGenerateUniqueIdentifier(mSeriesInstanceUID,SITE_SERIES_UID_ROOT); - dataset->putAndInsertString(DCM_SeriesInstanceUID, mSeriesInstanceUID); - dataset->putAndInsertString(DCM_FrameOfReferenceUID, mSeriesInstanceUID); - - ::initialDate(dataset, mStudyDate); - - dataset->putAndInsertString(DCM_StudyTime, mStudyTime.data()); - - ::initialSeriesTime(dataset); - - dataset->putAndInsertString(DCM_InstitutionName, mPatientData.getInstituationName().empty()? - "InstName":mPatientData.getInstituationName().data()); - dataset->putAndInsertString(DCM_InstitutionAddress, mPatientData.getInstituationAddress().empty()? - "default addr":mPatientData.getInstituationAddress().data()); - dataset->putAndInsertString(DCM_ReferringPhysicianName,mPatientData.getReferringPhysicianName().empty()? - "ReferringPhysician": mPatientData.getOperatorName().data()); - dataset->putAndInsertString(DCM_PatientName, mPatientData.getPatientName().empty()? - "TestPatient":mPatientData.getPatientName().data()); - dataset->putAndInsertString(DCM_PatientSex, mPatientData.getPatientSex().data()); - dataset->putAndInsertString(DCM_PatientBirthDate, mPatientData.getPatientBirthDate().data()); - dataset->putAndInsertString(DCM_PatientID, mPatientData.getPatientID().empty()? - "TestPatID":mPatientData.getPatientID().data()); - dataset->putAndInsertString(DCM_Laterality, mPatientData.getLaterality().data()); - - dataset->putAndInsertString(DCM_StudyID, "0"); - dataset->putAndInsertString(DCM_SeriesNumber, to_string(type).data()); - dataset->putAndInsertString(DCM_SeriesDescription, ::getSeriesDescription(type)); - // dataset->putAndInsertString(DCM_ProtocolName, "?"); - - // initial spacing and position - float spacing[3]{.0,.0,.0}; - spacing[0] = getSpacing(startPoint,endPoint,imageXYZ,0); - spacing[1] = getSpacing(startPoint,endPoint,imageXYZ,1); - spacing[2] = getSpacing(startPoint,endPoint,imageXYZ,2); - - float originPosition[3] = { - endPoint[1]*(float)1000.0, - endPoint[2]*(float)1000.0, - endPoint[0]*(float)1000.0, - }; - float originLocation =-endPoint[2]*1000.0; - dataset->putAndInsertString(DCM_SliceThickness, to_string(spacing[2]).data()); - dataset->putAndInsertUint16(DCM_Rows, Rows); - dataset->putAndInsertUint16(DCM_Columns, Cols); - string spacingsBP = to_string(spacing[0])+"\\"+to_string(spacing[1]); - dataset->putAndInsertString(DCM_PixelSpacing, spacingsBP.data()); - - // initial data, wwwl , Slope,Intercept - float Slope = 1, Intercept = 0 ; - size_t size = Rows*Cols*Slices; - ushort* udata = new ushort[size]{0}; - float min = data[0]; - float max = data[0]; - for (size_t i = 0; i < size; i++) - { - if (type == 0) data[i] = 1000*data[i]; - if (min > data[i]) min = data[i]; - if (max < data[i]) max = data[i]; - } - - long windowCenter = min + (max-min)/2, windowWidth = max-min; - float dSlope = 1.0/(pow(2,16)/windowWidth); - float dIntercept = min; - for (size_t i = 0; i < size; i++) - { - udata[i] = (ushort)((data[i] - min)/dSlope); - } - if (type !=0){ - Slope = dSlope; - Intercept = dIntercept; - } - else{ - windowCenter = (windowCenter-min)/dSlope; - windowWidth = (windowWidth-min)/dSlope; - } - dataset->putAndInsertString(DCM_WindowCenter, to_string(windowCenter).data()); - dataset->putAndInsertString(DCM_WindowWidth, to_string(windowWidth).data()); - dataset->putAndInsertString(DCM_RescaleSlope, to_string(Slope).data()); - dataset->putAndInsertString(DCM_RescaleIntercept, to_string(Intercept).data()); - - size_t sliceSize = Rows*Cols; - for (size_t i = 0; i < Slices; i++) - { - dataset->putAndInsertString(DCM_AccessionNumber, to_string(i).data()); - dataset->putAndInsertString(DCM_InstanceNumber, to_string(i).data()); - string pos = to_string(originPosition[1] ) + "\\" + - to_string(originPosition[2]+ i * spacing[2]) + "\\" + - to_string(originPosition[0]); - dataset->putAndInsertString(DCM_ImagePositionPatient, pos.data()); - string loc = to_string(originLocation + i * spacing[2]); - dataset->putAndInsertString(DCM_SliceLocation, loc.data()); - - // set SOPInstanceUID - char SOP_UID[100] = {0}; - dcmGenerateUniqueIdentifier(SOP_UID, SITE_INSTANCE_UID_ROOT); - metaInfo->putAndInsertString(DCM_MediaStorageSOPInstanceUID, SOP_UID); - dataset->putAndInsertString(DCM_SOPInstanceUID, SOP_UID); - - dataset->putAndInsertUint16Array(DCM_PixelData, udata + sliceSize * i, - sliceSize); - string path2 = path+to_string(i)+".dcm"; - OFFilename file ; - file = path2.data(); - auto status =dcmFile.saveFile(file, EXS_LittleEndianExplicit); - if(status.bad()) - { - RECON_ERROR("Save dicom file {0} fail, beacuse of {1}\r\n",path2.data(),status.text()); - } - else{ - RECON_INFO("Save dicom file {0} \r\n",path2.data()); - } - } - } -} \ No newline at end of file diff --git a/src/common/DICOMExporter.h b/src/common/DICOMExporter.h deleted file mode 100644 index 709805d..0000000 --- a/src/common/DICOMExporter.h +++ /dev/null @@ -1,41 +0,0 @@ -#ifndef __DICOMEXPORT_H__ -#define __DICOMEXPORT_H__ -#include - -#include "Matrix.h" - -#include "Parser/Data/PatientData.h" -namespace Recon { - -class DICOMExporter -{ -public: - enum ImageType{ - REFL, - SOS, - ATT - }; - explicit DICOMExporter(const PatientData& aPatientData); - //仅仅用于测试!!! - DICOMExporter(); - DICOMExporter(DICOMExporter &&) = default; - DICOMExporter(const DICOMExporter &) = default; - DICOMExporter &operator=(DICOMExporter &&) = default; - DICOMExporter &operator=(const DICOMExporter &) = default; - ~DICOMExporter(); - void setExportBasePath(const std::string & path); - void exportDICOM(Aurora::Matrix aMatrix, ImageType type); -private: - std::string mBasePath; - std::string mStudyDate; - std::string mStudyTime; - - char mStudyInstanceUID[100]={0}; - char mSeriesInstanceUID[100]={0}; - PatientData mPatientData; - void exportDCM(std::string path, float * startPoint, float* endPoint, int* imageXYZ,float* data, int type); -}; - - -} -#endif // __DICOMEXPORT_H__ \ No newline at end of file diff --git a/src/common/dataBlockCreation/getAScanBlockPreprocessed.cpp b/src/common/dataBlockCreation/getAScanBlockPreprocessed.cpp index f8bfb18..f428f61 100644 --- a/src/common/dataBlockCreation/getAScanBlockPreprocessed.cpp +++ b/src/common/dataBlockCreation/getAScanBlockPreprocessed.cpp @@ -90,70 +90,4 @@ AscanBlockPreprocessed Recon::getAscanBlockPreprocessed(Parser* aParser, const A } //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"< 0) - { - result.ascanBlockPreprocessed = preprocessAscanBlockCuda(result.ascanBlockPreprocessed, aMeasInfo); - } - // else - // { - // result.ascanBlockPreprocessed = ascanBlock.ascanBlock; - // } -//printTime(); -//std::cout<<"end"<getMetaData().getDataType() == "float16") { - ce = convertfp16tofloat(ceMeasuredPointer.get(), ceMeasuredLength, ceMeasuredPointer.getSize()); + ce = convertfp16tofloat(ceMeasuredPointer.get(), ceMeasuredLength, ceMeasuredPointer.getSize()); } else { diff --git a/src/common/preprocessAscanBlock.cpp b/src/common/preprocessAscanBlock.cpp index aad0f5e..5f7972e 100644 --- a/src/common/preprocessAscanBlock.cpp +++ b/src/common/preprocessAscanBlock.cpp @@ -21,14 +21,4 @@ Aurora::Matrix Recon::preprocessAscanBlock(const Aurora::Matrix& aAscans, const // end return result; -} - -Aurora::CudaMatrix Recon::preprocessAscanBlockCuda(const Aurora::CudaMatrix& aAscans, const MeasurementInfo& aMeasInfo) -{ - if(aMeasInfo.ascanDataType == "float16") - { - return Aurora::convertfp16tofloatCuda(aAscans, aAscans.getDimSize(0), aAscans.getDimSize(1)); - } - - return aAscans; } \ No newline at end of file diff --git a/src/common/preprocessAscanBlock.h b/src/common/preprocessAscanBlock.h index 30ed66c..b3563ec 100644 --- a/src/common/preprocessAscanBlock.h +++ b/src/common/preprocessAscanBlock.h @@ -8,7 +8,6 @@ namespace Recon { Aurora::Matrix preprocessAscanBlock(const Aurora::Matrix& aAscans, const MeasurementInfo& aMeasInfo); - Aurora::CudaMatrix preprocessAscanBlockCuda(const Aurora::CudaMatrix& aAscans, const MeasurementInfo& aMeasInfo); } #endif \ No newline at end of file diff --git a/src/config/config.cpp b/src/config/config.cpp index dd89991..5afba3b 100644 --- a/src/config/config.cpp +++ b/src/config/config.cpp @@ -307,11 +307,11 @@ namespace Recon nlohmann::json transmissionCorrection = params.at("transmissionCorrection"); if(transmissionCorrection.contains("soundSpeedCorrection")) { - reflectParams::soundSpeedCorrection = transmissionCorrection.at("soundSpeedCorrection").get(); + reflectParams::soundSpeedCorrection = 0; } if(transmissionCorrection.contains("attenuationCorrection")) { - reflectParams::attenuationCorrection = transmissionCorrection.at("attenuationCorrection").get(); + reflectParams::attenuationCorrection = 0; } if(transmissionCorrection.contains("saveTransmInReflCoords")) { @@ -679,7 +679,7 @@ namespace Recon reflectParams::windowLength = 10; reflectParams::numThreads = 30; reflectParams::expectedUSSpeedRange = Aurora::Matrix::fromRawData(new float[2] {1420,1600}, 1, 2); - //reflectParams.transmissionCorrection + // reflectParams.transmissionCorrection reflectParams::soundSpeedCorrection = 0; reflectParams::attenuationCorrection = 0; reflectParams::saveTransmInReflCoords = 1; diff --git a/src/main.cxx b/src/main.cxx index 73a261e..8b6e15b 100644 --- a/src/main.cxx +++ b/src/main.cxx @@ -12,7 +12,6 @@ #include "startReconstructions.h" #include "transmissionReconstruction/detection/detection.h" -#include "transmissionReconstruction/detection/detection.cuh" #include "startReconstructions.h" #include "log/log.h" /* 0 is data path. diff --git a/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp b/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp index bbc49a8..da64806 100644 --- a/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp +++ b/src/reflectionReconstruction/preprocessData/preprocessAScanBlockForReflection.cpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include namespace Recon { @@ -132,9 +133,11 @@ namespace Recon { } int winLength = 100; auto ascanMapValue = Aurora::zeros(1,blockedAScans.getDimSize(1),1); - #pragma omp parallel for num_threads(32) + // #pragma omp parallel for num_threads(32) + #pragma omp parallel for for (size_t i = 0; i < blockedAScans.getDimSize(1); i++) { + // mkl_set_num_threads_local(1); auto _blockedAScans = getPartMatrixColRef(blockedAScans,i); _blockedAScans = _blockedAScans/blockedGain[i]; if (reflectParams::removeDCOffset == 1) @@ -180,7 +183,8 @@ namespace Recon { if (reflectParams::useCorrelation && reflectParams::matchedFilterCeAScan) { - #pragma omp parallel for num_threads(32) + // #pragma omp parallel for num_threads(32) + #pragma omp parallel for for (size_t i = 0; i < blockedAScans.getDimSize(1); i++) { Aurora::Matrix fx = getPartMatrixColRef(fx_all,i); @@ -214,6 +218,7 @@ namespace Recon { Aurora::free(value2); fx_all(Aurora::$,i) = complex; } + blockedAScans = Aurora::real(ifft(fx_all)); blockedAScans.setBlockValue(0, 0, winLength-1, 0); @@ -225,15 +230,19 @@ namespace Recon { } auto exponent = offsetSignalFourier(preComputes, nSamples); // exponent.forceReshape(1, exponent.getDataSize(), 1); - #pragma omp parallel for num_threads(32) + + // #pragma omp parallel for num_threads(32) + #pragma omp parallel for for (size_t i = 0; i < blockedAScans.getDimSize(1); i++) { Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i); blockedAScans(Aurora::$,i) =_blockedAScans*exponent; } + blockedAScans = real(ifft(blockedAScans)); if (reflectParams::removeTransmissionSignal == 1) { - #pragma omp parallel for num_threads(32) + // #pragma omp parallel for num_threads(32) + #pragma omp parallel for for (size_t i = 0; i < blockedAScans.getDimSize(1); i++) { Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i); @@ -339,7 +348,8 @@ namespace Recon { if(reflectParams::useOptPulse==1){ auto n = nSamples; auto nHalf = round((float)n/2); - #pragma omp parallel for num_threads(32) + // #pragma omp parallel for num_threads(32) + #pragma omp parallel for for (size_t i = 0; i < blockedAScans.getDimSize(1); i++) { Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i); @@ -354,7 +364,9 @@ namespace Recon { blockedAScans = abs(blockedAScans); auto help_all = blockedAScans.deepCopy(); - #pragma omp parallel for num_threads(32) + + // #pragma omp parallel for num_threads(32) + #pragma omp parallel for for (size_t i = 0; i < blockedAScans.getDimSize(1); i++) { Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i); @@ -398,22 +410,26 @@ namespace Recon { help_all = fft(help_all); - #pragma omp parallel for num_threads(32) + // #pragma omp parallel for num_threads(32) + #pragma omp parallel for for (size_t i = 0; i < help_all.getDimSize(1); i++) { + // mkl_set_num_threads_local(1); Aurora::Matrix help = getPartMatrixColRef(help_all,i); help_all(Aurora::$,i) = help*(preComputes.sincPeak_ft); } help_all = real(ifft(help_all)); - #pragma omp parallel for num_threads(32) + // #pragma omp parallel for num_threads(32) + #pragma omp parallel for for (size_t i = 0; i < help_all.getDimSize(1); i++) { Aurora::Matrix _blockedAScans = getPartMatrixColRef(blockedAScans,i); Aurora::Matrix help = getPartMatrixColRef(help_all,i); Aurora::compareSet(_blockedAScans,0,help,Aurora::NL); } + } else{ blockedAScans = real(ifft(blockedAScans)); diff --git a/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp b/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp index 17da3fb..f35165e 100644 --- a/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp +++ b/src/reflectionReconstruction/reconstructionSAFT/reconstructionSAFT.cpp @@ -4,9 +4,6 @@ #include "Function2D.h" #include "Function3D.h" #include "Matrix.h" -#include "SAFTStructs.h" -// #include "SAFT_ATT.h" -#include "SAFT_TOFI.h" #include "common/dataBlockCreation/removeDataFromArrays.h" #include "config/config.h" @@ -19,314 +16,4 @@ #include namespace Recon { - - int compareColumn(const Aurora::Matrix &aMatrix, int colCount, - const Aurora::Matrix &aOther, int col2) { - - for (size_t j = 0; j < colCount; j++) { - size_t i = 0; - for (; i < aMatrix.getDimSize(0); i++) { - if (aMatrix[j * aMatrix.getDimSize(0) + 1] != - aOther[col2 * aOther.getDimSize(0) + i]) - break; - } - if (i == aMatrix.getDimSize(0)) - return j; - } - return -1; - } - - - Aurora::Matrix callSAFT(const Aurora::Matrix &AScans, - const Aurora::Matrix &receiverIndex, - const Aurora::Matrix &senderIndex, - const Aurora::Matrix &receiverPosGeom, - const Aurora::Matrix &senderPosGeom, - int SAFT_mode, float TimeInterval, - const TransRecos& transRecos, - const Aurora::Matrix &Env) - { - if (reflectParams::useAscanIndex == 0) { - if (reflectParams::soundSpeedCorrection == 0) { - std::cerr - << "reconstruction without sound speed correction not allowed!!!!" - << std::endl; - return Aurora::Matrix(); - } - } - std::vector params; - Matrix_t AScan{nullptr, - (size_t)AScans.getDims(), - {(size_t)AScans.getDimSize(0), (size_t)AScans.getDimSize(1), - (size_t)AScans.getDimSize(2)}, - AScans.getDataSize()}; - float *ASdata = new float[AScans.getDataSize()]{0}; - std::copy(AScans.getData(), AScans.getData() + AScans.getDataSize(), - ASdata); - AScan.Data = ASdata; - params.push_back(AScan); - - Matrix_t imageStartpoint{nullptr, 2, {1, 3}, 3}; - imageStartpoint.Data = - new float[3]{(float)reflectParams::imageStartpoint[0], - (float)reflectParams::imageStartpoint[1], - (float)reflectParams::imageStartpoint[2]}; - params.push_back(imageStartpoint); - - Matrix_t receiverIdx{nullptr, - 2, - { 1,receiverIndex.getDataSize(), 1}, - receiverIndex.getDataSize()}; - unsigned short *temp1 = new unsigned short[receiverIndex.getDataSize()]; - std::copy(receiverIndex.getData(), - receiverIndex.getData() + receiverIndex.getDataSize(), temp1); - receiverIdx.Data = temp1; - params.push_back(receiverIdx); - - Matrix_t senderIdx{nullptr, - 2, - {1,senderIndex.getDataSize(), 1}, - senderIndex.getDataSize()}; - auto temp2 = new unsigned short[senderIndex.getDataSize()]; - std::copy(senderIndex.getData(), - senderIndex.getData() + senderIndex.getDataSize(), temp2); - senderIdx.Data = temp2; - params.push_back(senderIdx); - - Matrix_t receiverGeom{nullptr, - (size_t)receiverPosGeom.getDims(), - {(size_t)receiverPosGeom.getDimSize(0), - (size_t)receiverPosGeom.getDimSize(1), - (size_t)receiverPosGeom.getDimSize(2)}, - receiverPosGeom.getDataSize()}; - float *fdata = new float[receiverPosGeom.getDataSize()]{0}; - std::copy(receiverPosGeom.getData(), - receiverPosGeom.getData() + receiverPosGeom.getDataSize(), fdata); - receiverGeom.Data = fdata; - params.push_back(receiverGeom); - - Matrix_t senderGeom{nullptr, - (size_t)senderPosGeom.getDims(), - {(size_t)senderPosGeom.getDimSize(0), - (size_t)senderPosGeom.getDimSize(1), - (size_t)senderPosGeom.getDimSize(2)}, - senderPosGeom.getDataSize()}; - auto fdata1 = new float[senderPosGeom.getDataSize()]{0}; - std::copy(senderPosGeom.getData(), - senderPosGeom.getData() + senderPosGeom.getDataSize(), fdata1); - senderGeom.Data = fdata1; - params.push_back(senderGeom); - - Matrix_t _SAFT_mode{&SAFT_mode, 1, {1, 1, 1}, 1}; - params.push_back(_SAFT_mode); - - int saftVariant_v[6]{(int)reflectParams::saftVariant[0],(int)reflectParams::saftVariant[1],(int)reflectParams::saftVariant[2], - (int)reflectParams::saftVariant[3],(int)reflectParams::saftVariant[4],(int)reflectParams::saftVariant[5]}; - Matrix_t saftVariant{saftVariant_v,2,{1,6,1},6}; - params.push_back(saftVariant); - - Matrix_t SpeedMap3D{nullptr, - (size_t)transRecos.speedMap3D.getDims(), - {(size_t)transRecos.speedMap3D.getDimSize(0), - (size_t)transRecos.speedMap3D.getDimSize(1), - (size_t)transRecos.speedMap3D.getDimSize(2)}, - transRecos.speedMap3D.getDataSize()}; - auto fdata2 = new float[transRecos.speedMap3D.getDataSize()]{0}; - std::copy(transRecos.speedMap3D.getData(), - transRecos.speedMap3D.getData() + transRecos.speedMap3D.getDataSize(), fdata2); - SpeedMap3D.Data = fdata2; - params.push_back(SpeedMap3D); - - Matrix_t BeginTransMap{nullptr, - (size_t)transRecos.beginTransMap.getDims(), - {(size_t)transRecos.beginTransMap.getDimSize(0), - (size_t)transRecos.beginTransMap.getDimSize(1), - (size_t)transRecos.beginTransMap.getDimSize(2)}, - transRecos.beginTransMap.getDataSize()}; - auto fdata3 = new float[transRecos.beginTransMap.getDataSize()]{0}; - std::copy(transRecos.beginTransMap.getData(), - transRecos.beginTransMap.getData() + transRecos.beginTransMap.getDataSize(), fdata3); - BeginTransMap.Data = fdata3; - params.push_back(BeginTransMap); - - Matrix_t DeltaTransMap{nullptr, - 1, - 1, - 1, - 1, - 1}; - auto fdata4 = new float[1]{0}; - fdata4[0] = transRecos.deltaTransMap; - DeltaTransMap.Data = fdata4; - params.push_back(DeltaTransMap); - - Matrix_t AttenuationMap3D{nullptr, - (size_t)transRecos.attenuationMap3D.getDims(), - {(size_t)transRecos.attenuationMap3D.getDimSize(0), - (size_t)transRecos.attenuationMap3D.getDimSize(1), - (size_t)transRecos.attenuationMap3D.getDimSize(2)}, - transRecos.attenuationMap3D.getDataSize()}; - auto fdata5 = new float[transRecos.attenuationMap3D.getDataSize()]{0}; - std::copy(transRecos.attenuationMap3D.getData(), - transRecos.attenuationMap3D.getData() + transRecos.attenuationMap3D.getDataSize(), fdata5); - AttenuationMap3D.Data = fdata5; - params.push_back(AttenuationMap3D); - - float imageResolution_v = (float)reflectParams::imageResolution; - Matrix_t imageResolution{&imageResolution_v, 1, {1, 1, 1}, 1}; - params.push_back(imageResolution); - - float TimeInterval_v = (float)TimeInterval; - Matrix_t _TimeInterval{&TimeInterval_v, 1, {1, 1, 1}, 1}; - params.push_back(_TimeInterval); - - int imageXYZ_v[3]{(int)reflectParams::imageXYZ[0],(int)reflectParams::imageXYZ[1],(int)reflectParams::imageXYZ[2]}; - Matrix_t imageXYZ{imageXYZ_v,2,{1,3,1},3}; - params.push_back(imageXYZ); - - Matrix_t _Env{nullptr, - (size_t)Env.getDims(), - {(size_t)Env.getDimSize(0), (size_t)Env.getDimSize(1), - (size_t)Env.getDimSize(2)}, - Env.getDataSize()}; - _Env.Data = Env.getData(); - params.push_back(_Env); - - int blockXYZ_v[3]{(int)reflectParams::blockDimXYZ[0],(int)reflectParams::blockDimXYZ[1],(int)reflectParams::blockDimXYZ[2]}; - Matrix_t blockXYZ{blockXYZ_v,2,{1,3,1},3}; - params.push_back(blockXYZ); - - Matrix_t gpus{nullptr, - 2, - {1, (size_t)reconParams::gpuSelectionList.getDataSize(), - 1}, - reconParams::gpuSelectionList.getDataSize()}; - int *gdata = new int[reconParams::gpuSelectionList.getDataSize()]{0}; - std::copy(reconParams::gpuSelectionList.getData(), reconParams::gpuSelectionList.getData() + reconParams::gpuSelectionList.getDataSize(), - gdata); - gpus.Data = gdata; - params.push_back(gpus); - - if (reflectParams::useAscanIndex == 0){ - Matrix_t medianWindowSize{&reflectParams::medianWindowSize,1,{1,1,1},1}; - params.push_back(medianWindowSize); - } - - float other[2]={(float)reflectParams::debugMode, (float)reflectParams::attenuationCorrectionLimit} ; - Matrix_t otherP{other,2,{1,2,1},2}; - params.push_back(otherP); - if (reflectParams::useAscanIndex == 0) { - // auto result = SAFT_ATT(params); - // delete [] ASdata; - // delete [] (float*)imageStartpoint.Data; - // delete [] temp1; - // delete [] temp2; - // delete [] fdata; - // delete [] fdata1; - // delete [] fdata2; - // delete [] fdata3; - // delete [] fdata4; - // delete [] fdata5; - // delete [] gdata; - // return Aurora::Matrix::fromRawData((float *)result.Data, result.Dims[0], - // result.Dims[1], result.Dims[2]); - - return Aurora::Matrix(); - } else { - auto result = SAFT_TOFI(params); - delete [] ASdata; - delete [] (float*)imageStartpoint.Data; - delete [] temp1; - delete [] temp2; - delete [] fdata; - delete [] fdata1; - delete [] fdata2; - delete [] fdata3; - delete [] fdata4; - delete [] fdata5; - delete [] gdata; - return Aurora::Matrix::fromRawData((float *)result.Data, result.Dims[0], - result.Dims[1], result.Dims[2]); - } - } - - Aurora::Matrix - 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) { - float TimeInterval = 1.0 / (reflectParams::aScanReconstructionFrequency); - std::vector motorPosAvailable; - std::unique_copy(mpIndex.getData(), mpIndex.getData() + mpIndex.getDataSize(), - std::back_inserter(motorPosAvailable)); - for (auto mp : motorPosAvailable) { - auto mpIdxs = mpIndex == mp; - size_t countMp = Aurora::sum(mpIdxs).getScalar(); - - auto senderPos_s = Aurora::transpose(removeDataFromArrays(senderPos, mpIdxs)); - auto receiverPos_s = Aurora::transpose(removeDataFromArrays(receiverPos, mpIdxs)); - Aurora::Matrix idsSender ; - auto senderPosGeom = Aurora::transpose(Aurora::uniqueByRows(senderPos_s, idsSender)); - idsSender = Aurora::transpose(idsSender); - Aurora::Matrix idsReceiver ; - auto receiverPosGeom = Aurora::transpose(Aurora::uniqueByRows(receiverPos_s, idsReceiver)); - idsReceiver = Aurora::transpose(idsReceiver); - - auto numUsedData = countMp; - int count = numUsedData / reflectParams::blockSizeReco; - auto blockIdxs = Aurora::zeros(1, count + 2); - for (size_t i = 0; i < count + 1; i++) { - blockIdxs[i] = i * reflectParams::blockSizeReco; - } - blockIdxs[count + 1] = numUsedData; - for (size_t i = 0; i < count + 1; i++) { - size_t length = blockIdxs[i + 1] - blockIdxs[i]; - size_t begin = blockIdxs[i]; - size_t end = blockIdxs[i + 1]; - - auto receiverIndex = Aurora::zeros(1, length); - auto senderIndex = Aurora::zeros(1, length); - auto mpIdxs_d = Aurora::zeros(1,mpIdxs.getDataSize()); - for (size_t j = begin, k = 0; j < end; j++, k++) { - mpIdxs_d[j] =mpIdxs[j]; - receiverIndex[k] = idsReceiver[j]+1; - senderIndex[k] = idsSender[j]+1; - } - auto _AScans = removeDataFromArrays(AScans, mpIdxs_d); - Env = callSAFT(_AScans, receiverIndex, senderIndex, receiverPosGeom, - senderPosGeom, SAFT_mode, TimeInterval, transRecos, Env); - } - } - return Env; - } - - //TODO:untested - Aurora::Matrix polyvaln(PolyModel polymodel, const Aurora::Matrix &indepvar) - { - auto np = Aurora::size(indepvar); - int n = (int)np[0]; - int p = (int)np[1]; - Aurora::Matrix _indepvar = indepvar; - if (n == 1 && polymodel.ModelTerms.getDimSize(1) == 1){ - _indepvar = Aurora::transpose(indepvar); - np = Aurora::size(indepvar); - n = (int)np[0]; - p = (int)np[1]; - } - else if(polymodel.ModelTerms.getDimSize(1) != p){ - std::cerr<<"Size of indepvar array and this model are inconsistent."< #include +#include #include #include @@ -39,18 +38,15 @@ void producerThread( Parser* aParser, const Aurora::Matrix& aMotorPos, const Aurora::Matrix& aRlList, const Aurora::Matrix& aRnList, GeometryInfo aGeom, MeasurementInfo aExpInfo, PreComputes aPreComputes) { -printf("Reflection reconstruction is carried out."); + // printf("Reflection reconstruction is carried out."); - printf("Preperations for reconstructions."); + // printf("Preperations for reconstructions."); - printf(" - reset GPUs"); + // printf(" - reset GPUs"); for (size_t i = 0; i < reflectParams::gpuSelectionList.getDataSize(); i++) { std::string msg; - if (!resetGPUDevice((int)reflectParams::gpuSelectionList[i],msg)) - { - std::cerr< lock(PRODUCER_MUTEX); - PRODUCER_BLOCKDATAS.push(blockData); - PRODUCER_PROCESSDATAS.push(preprocessData); - lock.unlock(); - PRODUCER_CONDITION.notify_one(); + std::cout<<"\rReflect reconstructing...... "<< index<<"/64"; + std::cout.flush(); + index++; + // std::unique_lock lock(PRODUCER_MUTEX); + // PRODUCER_BLOCKDATAS.push(blockData); + // PRODUCER_PROCESSDATAS.push(preprocessData); + // lock.unlock(); + // PRODUCER_CONDITION.notify_one(); } } } @@ -108,41 +113,42 @@ Aurora::Matrix Recon::startReflectionReconstruction( Parser* aParser, int aSAFT_ for (size_t i = 0; i < reflectParams::gpuSelectionList.getDataSize(); i++) { std::string msg; - if (!resetGPUDevice((int)reflectParams::gpuSelectionList[i],msg)) - { - std::cerr< lock(PRODUCER_MUTEX); + // PRODUCER_CONDITION.wait(lock, []{return !PRODUCER_PROCESSDATAS.empty() && !PRODUCER_BLOCKDATAS.empty();}); + // lock.unlock(); + // Env = recontructSAFT(removeDataFromArrays(PRODUCER_PROCESSDATAS.front().AscanBlock, PRODUCER_PROCESSDATAS.front().usedData), + // removeDataFromArrays(PRODUCER_BLOCKDATAS.front().senderPositionBlock, PRODUCER_PROCESSDATAS.front().usedData), + // removeDataFromArrays(PRODUCER_BLOCKDATAS.front().receiverPositionBlock, PRODUCER_PROCESSDATAS.front().usedData), + // removeDataFromArrays(PRODUCER_BLOCKDATAS.front().mpBlock, PRODUCER_PROCESSDATAS.front().usedData), + // aSAFT_mode, aTransRecos, Env); + // lock.lock(); + // PRODUCER_PROCESSDATAS.pop(); + // PRODUCER_BLOCKDATAS.pop(); + // lock.unlock(); - Matrix Env = Aurora::zeros((int)reflectParams::imageXYZ[0],(int)reflectParams::imageXYZ[1],(int)reflectParams::imageXYZ[2]); + // } + // Recon::notifyProgress(25+73*((j*i)/(aMotorPos.getDataSize() * aSlList.getDataSize()))); + // } + // } - for(int i=0; i lock(PRODUCER_MUTEX); - PRODUCER_CONDITION.wait(lock, []{return !PRODUCER_PROCESSDATAS.empty() && !PRODUCER_BLOCKDATAS.empty();}); - lock.unlock(); - Env = recontructSAFT(removeDataFromArrays(PRODUCER_PROCESSDATAS.front().AscanBlock, PRODUCER_PROCESSDATAS.front().usedData), - removeDataFromArrays(PRODUCER_BLOCKDATAS.front().senderPositionBlock, PRODUCER_PROCESSDATAS.front().usedData), - removeDataFromArrays(PRODUCER_BLOCKDATAS.front().receiverPositionBlock, PRODUCER_PROCESSDATAS.front().usedData), - removeDataFromArrays(PRODUCER_BLOCKDATAS.front().mpBlock, PRODUCER_PROCESSDATAS.front().usedData), - aSAFT_mode, aTransRecos, Env); - lock.lock(); - PRODUCER_PROCESSDATAS.pop(); - PRODUCER_BLOCKDATAS.pop(); - lock.unlock(); - std::cout< #include #include #include "Parser/Parser.h" @@ -33,11 +33,12 @@ using namespace Aurora; int Recon::startReconstructions(const std::string& aDataPath, const std::string& aDataRefPath, const std::string& aOutputPath) { - MatlabWriter writer(aOutputPath); + // omp_set_num_threads(32); + // mkl_set_num_threads(64); + + // MatlabWriter writer(aOutputPath); Parser dataParser(aDataPath); Parser refParser(aDataRefPath); - Recon::DICOMExporter exporter(dataParser.getPatientData()); - exporter.setExportBasePath(aOutputPath); if(!dataParser.getShotList()->isValid()) { RECON_INFO("startReconstructions with {0}.",aDataPath); @@ -224,8 +225,6 @@ int Recon::startReconstructions(const std::string& aDataPath, const std::string& 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"); } @@ -255,8 +254,7 @@ int Recon::startReconstructions(const std::string& aDataPath, const std::string& //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"); - exporter.exportDICOM(env, Recon::DICOMExporter::REFL); + // writer.setMatrix(env, "reflect"); // Recon::notifyProgress(99); } // writer.write(); diff --git a/src/transmissionReconstruction/detection/detection.cpp b/src/transmissionReconstruction/detection/detection.cpp index 22cff40..68f6c11 100644 --- a/src/transmissionReconstruction/detection/detection.cpp +++ b/src/transmissionReconstruction/detection/detection.cpp @@ -12,8 +12,6 @@ #include "common/getMeasurementMetaData.h" #include "config/config.h" -#include "calculateBankDetectAndHilbertTransformation.hpp" -#include "transmissionReconstruction/detection/detection.cuh" using namespace Aurora; namespace Recon { @@ -367,9 +365,6 @@ namespace Recon { resDetect = Aurora::zeros(1,N); resEnvelopeD = Aurora::malloc(totalSize); resEnvelopeRefD = Aurora::malloc(totalSize); - calculateBankDetectAndHilbertTransformation( - _AscanBlock_trim.getData(), _AscanRefBlock_trim.getData(), N, M, resampleFactor, nthreads, - resDetect.getData(), resEnvelopeD, resEnvelopeRefD); auto resEnvelope =Matrix::New(resEnvelopeD,M,N); auto resEnvelopeRef =Matrix::New(resEnvelopeRefD,M,N); //floor(size(AscanBlock, 1)*inits.resampleFactor / 2 - 1), @@ -415,15 +410,15 @@ 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( @@ -444,9 +439,9 @@ namespace Recon { 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(); + ret.att = r.att; + ret.sosValue = r.sosValue; + ret.tof = r.tof; return ret; } } diff --git a/src/transmissionReconstruction/detection/detection.cuh b/src/transmissionReconstruction/detection/detection.cuh_ similarity index 100% rename from src/transmissionReconstruction/detection/detection.cuh rename to src/transmissionReconstruction/detection/detection.cuh_ diff --git a/src/transmissionReconstruction/detection/detection.h b/src/transmissionReconstruction/detection/detection.h index 76a5292..b30a360 100644 --- a/src/transmissionReconstruction/detection/detection.h +++ b/src/transmissionReconstruction/detection/detection.h @@ -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 diff --git a/src/transmissionReconstruction/detection/getTransmissionData.cpp b/src/transmissionReconstruction/detection/getTransmissionData.cpp index b78718e..295ce74 100644 --- a/src/transmissionReconstruction/detection/getTransmissionData.cpp +++ b/src/transmissionReconstruction/detection/getTransmissionData.cpp @@ -1,11 +1,9 @@ #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/log.h" #include "log/notify.h" #include "src//config/config.h" #include "src/common/getGeometryInfo.h" @@ -31,11 +29,6 @@ #include #include -#include -#include "Function1D.cuh" -#include "Function2D.cuh" -#include - using namespace Recon; using namespace Aurora; @@ -50,8 +43,8 @@ namespace Matrix waterTempBlock; MetaInfos metaInfos; - CudaMatrix ascanBlock; - CudaMatrix ascanBlockRef; + Matrix ascanBlock; + Matrix ascanBlockRef; Matrix dists; Matrix distRefBlock; Matrix waterTempRefBlock; @@ -64,59 +57,28 @@ namespace std::mutex PROCESS_BUFFER_MUTEX; std::condition_variable PROCESS_BUFFER_CONDITION; int BUFFER_COUNT = 0; - int BUFFER_SIZE = 4;//<=8 + int BUFFER_SIZE = 4; -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 << ":" < lock(CREATE_BUFFER_MUTEX); BLOCK_OF_TRANSIMISSIONDARA_BUFFER[std::to_string(aIndex)] = buffer; - std::cout<<"Add: "< lock(PROCESS_BUFFER_MUTEX); @@ -354,9 +349,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,24 +370,28 @@ 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 lockBufferCount(CREATE_BUFFER_MUTEX); BLOCK_OF_TRANSIMISSIONDARA_BUFFER.erase(std::to_string(index)); --BUFFER_COUNT; + finish_count++; lockBufferCount.unlock(); - std::cout<<"Remove: "< potentialMapDataDims = {1,1,1}; for(size_t i=0; i namespace Recon { - bool isEigsFinished = false; - bool isEigsStatus = false; - - float eigs(Aurora::Sparse& aM) - { - size_t size = aM.getM(); - if(size < aM.getN()) - { - size = aM.getN(); - } - Eigen::SparseMatrix M(size,size); - Aurora::Matrix rows = aM.getRowVector(); - Aurora::Matrix columns = aM.getColVector(); - Aurora::Matrix values = aM.getValVector(); - std::vector> triplets(rows.getDataSize()); - #pragma omp parallel for - for (int i = 0; i < rows.getDataSize(); ++i) - { - triplets[i] = Eigen::Triplet(rows[i], columns[i], values[i]); - } - M.setFromTriplets(triplets.begin(), triplets.end()); - float result; - Spectra::SparseGenMatProd op(M); - Spectra::GenEigsSolver> eigs(op, 1, 6); - eigs.init(); - int nconv = eigs.compute(Spectra::SortRule::LargestMagn); - Eigen::VectorXcd evalues; - if(eigs.info() == Spectra::CompInfo::Successful) - { - evalues = eigs.eigenvalues(); - std::complex complex = evalues[0]; - result = complex.real(); - } - isEigsFinished = true; - if (result> 1 + 1e-10) - { - isEigsStatus = true; - } - return result; - } - - void checkAndScale(Aurora::Sparse& M, Aurora::Matrix& b,size_t n) - { - float s2 = 0.0; - if(!isEigsFinished) - { - s2 = eigs(M); - return;; - } - - if (isEigsStatus) - { - b = b / sqrt(s2); - M.getValVector() = M.getValVector() / sqrt( s2); - } - } - - Aurora::Matrix callTval3(Aurora::Sparse& M, Aurora::Matrix& b,const Aurora::Matrix& dims,int device, TVALOptions& opt) - { - checkAndScale(M,b,(size_t)Aurora::prod(dims).getScalar()); - int * yIdxs = new int[M.getColVector().getDataSize()]; - std::copy(M.getColVector().getData(),M.getColVector().getData()+M.getColVector().getDataSize(),yIdxs); - int * xIdxs = new int[M.getRowVector().getDataSize()]; - std::copy(M.getRowVector().getData(),M.getRowVector().getData()+M.getRowVector().getDataSize(),xIdxs); - Aurora::Matrix values = M.getValVector(); - size_t cols = M.getM(), rows = M.getN(); - int nz = std::max(M.getColVector().getDataSize(),std::max(M.getRowVector().getDataSize(),M.getValVector().getDataSize())); - sparse_matrix_t A; - sparse_matrix_t csrA; - mkl_sparse_s_create_coo(&A, sparse_index_base_t::SPARSE_INDEX_BASE_ZERO, rows, cols, nz, yIdxs, xIdxs,values.getData()); - mkl_sparse_convert_csr(A, sparse_operation_t::SPARSE_OPERATION_NON_TRANSPOSE, &csrA); - int n_rows,n_cols; - int *rows_start,*rows_end,*col_indx; - float * csrValues; - sparse_index_base_t index; - mkl_sparse_s_export_csr(csrA, &index, &n_rows, &n_cols, &rows_start, &rows_end, &col_indx, &csrValues); - mkl_sparse_destroy(A); - delete [] xIdxs; - delete [] yIdxs; - - int *row_idx = new int[n_rows+1]; - std::copy(rows_start,rows_start+n_rows,row_idx); - row_idx[n_rows] = rows_end[n_rows-1]; - float* bData = b.getData(); - // std::copy(b.getData(),b.getData()+b.getDataSize(),bData); - size_t bDims[3]={(size_t)b.getDimSize(0),(size_t)b.getDimSize(1),(size_t)b.getDimSize(2)}; - size_t rdims[3] = {(size_t)dims[0], (size_t)dims[1], (size_t)dims[2]}; - bool pagelocked = false; - auto result = TVALGPU( col_indx, row_idx, csrValues, M.getM(), M.getN(), nz, bData, bDims, rdims, opt, device, false); - mkl_sparse_destroy(csrA); - delete [] row_idx; - //delete [] bData; - return Aurora::Matrix::fromRawData(result.data, result.dims[0],result.dims[1],result.dims[2]); - // return Aurora::Matrix(); - } } diff --git a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h index 64b634a..e4b7e27 100644 --- a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h +++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h @@ -1,17 +1,8 @@ #ifndef __TVAL_H__ #define __TVAL_H__ #include "Sparse.h" -#include "tvalstruct.h" namespace Recon { - -void checkAndScale(Aurora::Sparse& M, Aurora::Matrix& b,size_t n); - -float eigs(Aurora::Sparse& aM); - -Aurora::Matrix callTval3(Aurora::Sparse &M, Aurora::Matrix &b, - const Aurora::Matrix &dims, int device, - struct TVALOptions &options); } #endif // __TVAL_H__ \ No newline at end of file diff --git a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.cpp b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.cpp index d447210..e03141d 100644 --- a/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.cpp +++ b/src/transmissionReconstruction/reconstruction/solvingEquationSystem/solve.cpp @@ -8,7 +8,6 @@ #include "config/config.h" #include "transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h" -#include "tvalstruct.h" namespace Recon { @@ -26,25 +25,7 @@ namespace Recon if (Recon::transParams::name.empty()){ Recon::transParams::name = "TVAL3"; } - if (Recon::transParams::name == "TVAL3") { - //callTval3 - TVALOptions opt; - opt.nonneg = solverOptions.nonNeg; - opt.bent = false; - opt.tol = 1E-10; - opt.maxit = niter; - opt.TVnorm = 2; - opt.disp = false; - opt.mu0 = solverOptions.TVAL3MU0; - opt.mu = solverOptions.TVAL3MU; - opt.beta = solverOptions.TVAL3Beta; - opt.beta0 = solverOptions.TVAL3Beta0; - int device = aDevice; - return callTval3(M, b, dims, device, opt); - } - //SART - else{ //TODO:待实现,先实现默认的TVAL3 return Aurora::Matrix(); }