UR/src/common/DICOMExporter.cpp

#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  {
    float getSpacing(float * startPoint, float* endPoint, int* imageXYZ, int index){
    return (endPoint[index]-startPoint[index])/(float)(imageXYZ[index])*1000.0;
}

    void initialDicomFile(DcmDataset* aDataset,DcmMetaInfo* aMetaInfo)
    {
        aMetaInfo->putAndInsertString(DCM_MediaStorageSOPClassUID,
                                        UID_CTImageStorage);
        aDataset->putAndInsertString(DCM_SOPClassUID, UID_CTImageStorage);

        aDataset->putAndInsertString(DCM_TransferSyntaxUID,
                                        UID_LittleEndianImplicitTransferSyntax);
        aDataset->putAndInsertString(DCM_ImageType, "ORIGINAL\\PRIMARY");

        aDataset->putAndInsertString(DCM_Manufacturer, "EQ9");
        aDataset->putAndInsertString(DCM_PatientOrientation, "L\\H");
        aDataset->putAndInsertString(DCM_PatientPosition, "HFP");

        aDataset->putAndInsertString(DCM_ImageOrientationPatient, "1\\0\\0\\0\\0\\1");

        aDataset->putAndInsertString(DCM_StudyDescription, "USCT BREAST");
        aDataset->putAndInsertString(DCM_ManufacturerModelVersion, "USCT3DV3");

        aDataset->putAndInsertString(DCM_SamplesPerPixel, "1");
        aDataset->putAndInsertString(DCM_PixelRepresentation, "0");
        aDataset->putAndInsertString(DCM_BitsAllocated, "16");
        aDataset->putAndInsertString(DCM_BitsStored, "16");
        aDataset->putAndInsertString(DCM_HighBit, "15");
        aDataset->putAndInsertString(DCM_PhotometricInterpretation, "MONOCHROME2");

        aDataset->putAndInsertString(DCM_KVP, "0");

    }
    void initialDate(DcmDataset* aDataset,const std::string& aStudyDate )
    {
        aDataset->putAndInsertString(DCM_StudyDate, aStudyDate.data());
        OFDateTime currentDateTime;
        currentDateTime.setCurrentDateTime();
        OFString date;
        currentDateTime.getDate().getISOFormattedDate(date,false);
        aDataset->putAndInsertString(DCM_SeriesDate, aStudyDate.data());
        aDataset->putAndInsertString(DCM_AcquisitionDate, aStudyDate.data());
        aDataset->putAndInsertString(DCM_ContentDate, date.data());
    }

    void initialSeriesTime(DcmDataset* aDataset,const std::string& aStudyTime){
                // 获取当前日期和时间（UTC时间）
        OFDateTime currentDateTime;
        currentDateTime.setCurrentDateTime();
        OFString time;
        currentDateTime.getTime().getISOFormattedTime(time,true,false,false,false);
        aDataset->putAndInsertString(DCM_SeriesTime, aStudyTime.data());
        aDataset->putAndInsertString(DCM_AcquisitionTime, aStudyTime.data());
        aDataset->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, const MetaData& aMetaData)
    {
        mPatientData = aPatientData;
        //prepare Study
        if (aPatientData.getStudyUID().empty()){
            char tempStr [100] = {0};
            dcmGenerateUniqueIdentifier(tempStr,SITE_STUDY_UID_ROOT);
            mStudyInstanceUID.append(tempStr);
        }
        else{
            mStudyInstanceUID = aPatientData.getStudyUID();
        }
        std::string measurementID = aMetaData.getMeasurementID();
        
        unsigned long dateIndex = measurementID.find('_');
        unsigned long timeIndex = measurementID.find('T');
        OFDateTime currentDateTime;
        currentDateTime.setCurrentDateTime();
        if(dateIndex== std::string::npos || measurementID.length()-dateIndex<=8){

            OFString date;
            currentDateTime.getDate().getISOFormattedDate(date,false);
            mStudyDate = date.data();

        }
        else{
            mStudyDate = measurementID.substr(dateIndex+1,8);
        }
        if(timeIndex== std::string::npos || measurementID.length()-timeIndex<=6){
            OFString time;
            currentDateTime.getTime().getISOFormattedTime(time,true,false,false,false);
            mStudyTime = time.data();
        }
        else{
            mStudyTime = measurementID.substr(timeIndex+1,6);
        }
    }
    
    DICOMExporter::DICOMExporter()
    {
        //prepare Study
        char tempStr [100] = {0};
        dcmGenerateUniqueIdentifier(tempStr,SITE_STUDY_UID_ROOT);
        mStudyInstanceUID.append(tempStr);

        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_SpecificCharacterSet, "ISO_IR 192");

        dataset->putAndInsertString(DCM_StudyInstanceUID, mStudyInstanceUID.data());

        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, mStudyTime);

        dataset->putAndInsertString(DCM_Modality, mPatientData.getModality().empty()?
            "CT":mPatientData.getModality().data());

        dataset->putAndInsertString(DCM_InstitutionName, mPatientData.getInstituationName().empty()?
            "":mPatientData.getInstituationName().data());
        dataset->putAndInsertString(DCM_InstitutionAddress, mPatientData.getInstituationAddress().empty()?
            "":mPatientData.getInstituationAddress().data());
        dataset->putAndInsertString(DCM_ReferringPhysicianName,mPatientData.getReferringPhysicianName().empty()?
            "": mPatientData.getReferringPhysicianName().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, mStudyTime.data());
        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, mPatientData.getAccessionNumber().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());
            }
        }
    }
}