Add DICOMExporter and use it to save images.

src/common/DICOMExporter.cpp

@@ -0,0 +1,239 @@
+#include "DICOMExporter.h"
+
+#include "fileHelper.h"
+#include "config/config.h"
+#include "log/log.h"
+
+#include <cstddef>
+#include <cstdio>
+#include <dcmtk/dcmdata/dcdeftag.h>
+#include <dcmtk/dcmdata/dcfilefo.h>
+#include <dcmtk/dcmdata/dcuid.h>
+#include <dcmtk/ofstd/offile.h>
+#include <string>
+#include <sys/types.h>
+#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());
+            }
+        }
+    }
+}

src/common/DICOMExporter.h

@@ -0,0 +1,39 @@
+#ifndef __DICOMEXPORT_H__
+#define __DICOMEXPORT_H__
+#include <string>
+
+#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__

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");
+        exporter.exportDICOM(transmissionResult.recoSOS, DICOMExporter::SOS);
-        writer.setMatrix(transmissionResult.recoATT, "att");
+        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();
 }