diff --git a/src/common/DICOMExporter.cpp b/src/common/DICOMExporter.cpp new file mode 100644 index 0000000..73f8917 --- /dev/null +++ b/src/common/DICOMExporter.cpp @@ -0,0 +1,239 @@ +#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 { + double getSpacing(double * startPoint, double* endPoint, int* imageXYZ, int index){ + return (endPoint[index]-startPoint[index])/(double)(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() + { + } + 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); + } + + void DICOMExporter::exportDCM(string path, double * startPoint, double* endPoint, int* imageXYZ,double* data, int type) + { + long Rows = imageXYZ[0], Cols = imageXYZ[1] ,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().data()); + dataset->putAndInsertString(DCM_InstitutionAddress, mPatientData.getInstituationAddress().data()); + dataset->putAndInsertString(DCM_ReferringPhysicianName, mPatientData.getOperatorName().data()); + dataset->putAndInsertString(DCM_PatientName, mPatientData.getPatientName().data()); + dataset->putAndInsertString(DCM_PatientSex, mPatientData.getPatientSex().data()); + dataset->putAndInsertString(DCM_PatientBirthDate, mPatientData.getPatientBirthDate().data()); + dataset->putAndInsertString(DCM_PatientID, 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 + double 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); + + double originPosition[3] = { + endPoint[1]*1000.0, + endPoint[2]*1000.0, + endPoint[0]*1000.0, + }; + double originLocation =endPoint[1]*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 + double Slope = 1, Intercept = 0 ; + size_t size = Rows*Cols*Slices; + ushort* udata = new ushort[size]{0}; + double min = data[0]; + double 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; + double dSlope = 1.0/(pow(2,16)/windowWidth); + double 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[0] ) + "\\" + + to_string(originPosition[1]- i * spacing[2]) + "\\" + + to_string(originPosition[2]); + 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 new file mode 100644 index 0000000..82e1314 --- /dev/null +++ b/src/common/DICOMExporter.h @@ -0,0 +1,39 @@ +#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 &&) = 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, double * startPoint, double* endPoint, int* imageXYZ,double* data, int type); +}; + + +} +#endif // __DICOMEXPORT_H__ \ No newline at end of file diff --git a/src/startReconstructions.cpp b/src/startReconstructions.cpp index 58798ae..1325cbf 100644 --- a/src/startReconstructions.cpp +++ b/src/startReconstructions.cpp @@ -4,6 +4,7 @@ #include "MatlabWriter.h" #include "config/config.h" #include "log/log.h" +#include "common/DICOMExporter.h" #include "common/getMeasurementMetaData.h" #include "common/getGeometryInfo.h" #include "common/estimatePulseLength.h" @@ -34,6 +35,8 @@ void Recon::startReconstructions(const std::string& aDataPath, const std::string MatlabWriter writer(aOutputPath); Parser dataParser(aDataPath); Parser refParser(aDataRefPath); + Recon::DICOMExporter exporter(dataParser.getPatientData()); + exporter.setExportBasePath(aOutputPath); if(!dataParser.getShotList()->isValid()) { RECON_INFO("data path is invalid."); @@ -219,8 +222,10 @@ void 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; - writer.setMatrix(transmissionResult.recoSOS, "sos"); - writer.setMatrix(transmissionResult.recoATT, "att"); + exporter.exportDICOM(transmissionResult.recoSOS, DICOMExporter::SOS); + exporter.exportDICOM(transmissionResult.recoATT, DICOMExporter::ATT); + // writer.setMatrix(transmissionResult.recoSOS, "sos"); + // writer.setMatrix(transmissionResult.recoATT, "att"); } if(reflectParams::runReflectionReco) @@ -245,8 +250,9 @@ void Recon::startReconstructions(const std::string& aDataPath, const std::string reflectParams::gpuSelectionList = reconParams::gpuSelectionList; 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); } - writer.write(); + // writer.write(); } \ No newline at end of file