feat: Add position image in DICOMExporter.

src/common/DICOMExporter.cpp

@@ -11,6 +11,7 @@
 #include <dcmtk/dcmdata/dctk.h>
 
 #include "fileHelper.h"
+#include "Function1D.h"
 #include "config/config.h"
 #include "log/log.h"
 
@@ -20,35 +21,67 @@ namespace  {
     return (endPoint[index]-startPoint[index])/(float)(imageXYZ[index])*1000.0;
 }
 
-    void initialDicomFile(DcmDataset* aDataset,DcmMetaInfo* aMetaInfo)
+    Aurora::Matrix flipHorizontal(const Aurora::Matrix& aImage)
     {
-        aMetaInfo->putAndInsertString(DCM_MediaStorageSOPClassUID,
+        int row = aImage.getDimSize(0);
-                                        UID_CTImageStorage);
+        int column = aImage.getDimSize(1);
-        aDataset->putAndInsertString(DCM_SOPClassUID, UID_CTImageStorage);
+        Aurora::Matrix result = Aurora::Matrix::fromRawData(new float[row*column], row, column);
-
+        for (int c = 0; c < column; ++c)
-        aDataset->putAndInsertString(DCM_TransferSyntaxUID,
+        {
-                                        UID_LittleEndianImplicitTransferSyntax);
+            const int mirrorC = column - 1 - c;
-        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");
             
+            for (int r = 0; r < row; ++r)
+            {
+                const int originalIndex = c * row + r;
+                const int newIndex = mirrorC * row + r;
+                result[newIndex] = aImage[originalIndex];
             }
+        } 
+        return result;
+    }
+
+    Aurora::Matrix rotate90Clock2D(const Aurora::Matrix& aImage)
+    {
+        int row = aImage.getDimSize(0);
+        int column = aImage.getDimSize(1);
+        Aurora::Matrix result = Aurora::Matrix::fromRawData(new float[row*column], column, row);
+        const int originalRow = row;
+        const int originalCol = column;
+        
+        const int newRow = originalCol; 
+        const int newCol = originalRow;  
+        for (int newColIndex = 0; newColIndex < newCol; ++newColIndex)
+        {
+            for (int newRowIndex = 0; newRowIndex < newRow; ++newRowIndex)
+            {
+                const int oldColindex = newRowIndex;  
+                const int oldRowindex = originalRow - 1 - newColIndex;  
+                
+                const int oldPos = oldColindex * originalRow + oldRowindex;
+                const int newPos = newColIndex * newRow + newRowIndex;
+                
+                result[newPos] = aImage[oldPos];
+            }
+        }
+        return result;
+    }
+
+    Aurora::Matrix transferToDicomType(const Aurora::Matrix& aImage, Recon::DICOMExporter::ImageType aType)
+    {
+        switch (aType)
+        {
+        case Recon::DICOMExporter::REFL:
+        case Recon::DICOMExporter::CoronalPositioning:       
+            return Aurora::transpose(flipHorizontal(aImage));
+        case Recon::DICOMExporter::SagittalPositioning:
+            return Aurora::transpose(aImage);
+        case Recon::DICOMExporter::TransversePositioning:
+            return  Aurora::transpose(rotate90Clock2D(aImage)); 
+        default:
+            return aImage;
+        }
+    }
+
     void initialDate(DcmDataset* aDataset,const std::string& aStudyDate )
     {
         aDataset->putAndInsertString(DCM_StudyDate, aStudyDate.data());
@@ -71,7 +104,7 @@ namespace  {
         aDataset->putAndInsertString(DCM_AcquisitionTime, aStudyTime.data());
         aDataset->putAndInsertString(DCM_ContentTime, time.data());
     }
-    const char* SeriesDescriptions[3]={"REFL","SOS","ATT"};
+    const char* SeriesDescriptions[6]={"REFL","SOS","ATT","CoronalPositioning","SagittalPositioning","TransversePositioning"};
     const char* getSeriesDescription(int type){
         return SeriesDescriptions[type];
     }
@@ -83,6 +116,7 @@ namespace Recon
     DICOMExporter::DICOMExporter(const PatientData& aPatientData, const MetaData& aMetaData)
     {
         mPatientData = aPatientData;
+        mDicomFile = new DcmFileFormat();
         //prepare Study
         if (aPatientData.getStudyUID().empty()){
             char tempStr [100] = {0};
@@ -93,7 +127,7 @@ namespace Recon
             mStudyInstanceUID = aPatientData.getStudyUID();
         }
         std::string measurementID = aMetaData.getMeasurementID();
-        
+        mLaterality = mPatientData.getLaterality().empty() ? "L" : mPatientData.getLaterality();
         unsigned long dateIndex = measurementID.find('_');
         unsigned long timeIndex = measurementID.find('T');
         OFDateTime currentDateTime;
@@ -116,10 +150,12 @@ namespace Recon
         else{
             mStudyTime = measurementID.substr(timeIndex+1,6);
         }
+        initialDicomFile();
     }
     
     DICOMExporter::DICOMExporter()
     {
+        mDicomFile = new DcmFileFormat();
         //prepare Study
         char tempStr [100] = {0};
         dcmGenerateUniqueIdentifier(tempStr,SITE_STUDY_UID_ROOT);
@@ -133,44 +169,48 @@ namespace Recon
         currentDateTime.getTime().getISOFormattedTime(time,true,false,false,false);
         mStudyDate = date.data();
         mStudyTime = time.data();
+        initialDicomFile();
     }
 
     DICOMExporter::~DICOMExporter()
     {
+        delete mDicomFile;
     }
     void DICOMExporter::setExportBasePath(const std::string & path)
     {
         mBasePath = path;
     }
     
-    void DICOMExporter::exportDICOM(Aurora::Matrix aMatrix, ImageType type){
+    void DICOMExporter::initialDicomFile()
-        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];
+        DcmDataset* dataset = mDicomFile->getDataset();
-        DcmFileFormat dcmFile;
+        DcmMetaInfo* metaInfo =  mDicomFile->getMetaInfo();
-        DcmDataset* dataset = dcmFile.getDataset();
+        metaInfo->putAndInsertString(DCM_MediaStorageSOPClassUID,
-        DcmMetaInfo* metaInfo = dcmFile.getMetaInfo();
+                                        UID_CTImageStorage);
+        dataset->putAndInsertString(DCM_SOPClassUID, UID_CTImageStorage);
+        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_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");
 
-        initialDicomFile(dataset,metaInfo);
         dataset->putAndInsertString(DCM_SpecificCharacterSet, "ISO_IR 192");
 
         dataset->putAndInsertString(DCM_StudyInstanceUID, mStudyInstanceUID.data());
 
-        dcmGenerateUniqueIdentifier(mSeriesInstanceUID,SITE_SERIES_UID_ROOT);
+        dcmGenerateUniqueIdentifier(mFrameOfReferenceUID,SITE_SERIES_UID_ROOT);
-        dataset->putAndInsertString(DCM_SeriesInstanceUID, mSeriesInstanceUID);
+        dataset->putAndInsertString(DCM_FrameOfReferenceUID, mFrameOfReferenceUID);
-        dataset->putAndInsertString(DCM_FrameOfReferenceUID, mSeriesInstanceUID);
 
         ::initialDate(dataset, mStudyDate);
 
@@ -196,32 +236,148 @@ namespace Recon
         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, "?");
 
+    void DICOMExporter::initImageOrientation(ImageType aType)
+    {
+        DcmDataset* dataset = mDicomFile->getDataset();
+        DcmMetaInfo* metaInfo =  mDicomFile->getMetaInfo();
+        dataset->putAndInsertString(DCM_ImageOrientationPatient, getImageOrientationPatient(aType).c_str());
+        dataset->putAndInsertString(DCM_SeriesNumber, to_string(aType).data());
+        dataset->putAndInsertString(DCM_SeriesDescription, std::string(std::string(::getSeriesDescription(aType)) + "-" + mLaterality).c_str());
+    }
+
+    std::string DICOMExporter::getImageOrientationPatient(ImageType aType)
+    {
+        switch (aType)
+        {
+        case REFL:
+            return "1\\0\\0\\0\\0\\-1";
+        case SOS:
+            return "1\\0\\0\\0\\0\\-1";
+        case ATT:
+            return "1\\0\\0\\0\\0\\-1";
+        case CoronalPositioning:
+            return "1\\0\\0\\0\\0\\-1";
+        case SagittalPositioning:
+            return "0\\1\\0\\0\\0\\-1";
+        case TransversePositioning:
+            return "1\\0\\0\\0\\1\\0"; 
+        }
+    }
+
+    void DICOMExporter::initOriginPosition(ImageType aType, float* aOriginPosition, float* aPoint, int* aXYZ, float* aSpacing)
+    {
+        switch (aType)
+        {
+        case REFL:
+        case SOS:
+        case ATT:
+            aOriginPosition[0] = aPoint[1]*(float)1000.0;
+            aOriginPosition[1] = aPoint[2]*(float)1000.0;
+            aOriginPosition[2] = aPoint[0]*(float)1000.0;
+            break;
+        case CoronalPositioning:
+            aOriginPosition[0] = aPoint[1]*(float)1000.0;
+            aOriginPosition[1] = aPoint[2]*(float)1000.0 + (aXYZ[2] / 2 - 1) * aSpacing[2];
+            aOriginPosition[2] = aPoint[0]*(float)1000.0;
+            break;
+            case SagittalPositioning:
+            aOriginPosition[0] = aPoint[1]*(float)1000.0 + (aXYZ[1] / 2 - 1)* aSpacing[1];
+            aOriginPosition[1] = aPoint[2]*(float)1000.0;
+            aOriginPosition[2] = aPoint[0]*(float)1000.0;
+            break;
+            case TransversePositioning:
+            aOriginPosition[0] = aPoint[1]*(float)1000.0;
+            aOriginPosition[1] = aPoint[2]*(float)1000.0;
+            aOriginPosition[2] = aPoint[0]*(float)1000.0 - (aXYZ[0] / 2 - 1) * aSpacing[0];
+            break;
+        default:
+            break;
+        }
+    }
+
+    void DICOMExporter::initPixelSpacing(ImageType aType, float* aSpacing, int aRow, int aColumn)
+    {
+        DcmDataset* dataset = mDicomFile->getDataset();
+        DcmMetaInfo* metaInfo =  mDicomFile->getMetaInfo();
+        string spacingsBP;
+        switch (aType)
+        {
+        case REFL:
+        case SOS:
+        case ATT:
+        case CoronalPositioning:
+        dataset->putAndInsertString(DCM_SliceThickness, to_string(aSpacing[2]).data());
+        dataset->putAndInsertUint16(DCM_Rows, aRow);
+        dataset->putAndInsertUint16(DCM_Columns, aColumn);
+        spacingsBP = to_string(aSpacing[1])+"\\"+to_string(aSpacing[0]);
+        dataset->putAndInsertString(DCM_PixelSpacing, spacingsBP.data());
+        break;
+        case SagittalPositioning:
+        dataset->putAndInsertString(DCM_SliceThickness, to_string(aSpacing[1]).data());
+        dataset->putAndInsertUint16(DCM_Rows, aRow);
+        dataset->putAndInsertUint16(DCM_Columns, aColumn);
+        spacingsBP = to_string(aSpacing[2])+"\\"+to_string(aSpacing[0]);
+        dataset->putAndInsertString(DCM_PixelSpacing, spacingsBP.data());
+        break;
+        case TransversePositioning:
+        dataset->putAndInsertString(DCM_SliceThickness, to_string(aSpacing[0]).data());
+        dataset->putAndInsertUint16(DCM_Rows, aColumn);//因为旋转，必须交换row和column
+        dataset->putAndInsertUint16(DCM_Columns, aRow);
+        spacingsBP = to_string(aSpacing[1])+"\\"+to_string(aSpacing[2]);
+        dataset->putAndInsertString(DCM_PixelSpacing, spacingsBP.data());
+        break;
+        }
+    }
+
+    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, XYZ);
+        RECON_INFO("Save DICOM to Path:{0}, type:{1}==================>",path,(int)type);
+    }
+
+    void DICOMExporter::exportDICOM(Aurora::Matrix aMatrix, const Aurora::Matrix& aOriginMatrix, ImageType type)
+    {
+        int XYZ [3] = {aMatrix.getDimSize(0),aMatrix.getDimSize(1), aMatrix.getDimSize(2)};
+        int originXYZ[3] = {aOriginMatrix.getDimSize(0),aOriginMatrix.getDimSize(1), aOriginMatrix.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, originXYZ);
+        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, ImageType type, int* aOriginXYZ)
+    {
+        long Cols = imageXYZ[1], Rows = imageXYZ[0] ,Slices = imageXYZ[2];
+        // dataset->putAndInsertString(DCM_ProtocolName, "?");
+        DcmDataset* dataset = mDicomFile->getDataset();
+        DcmMetaInfo* metaInfo =  mDicomFile->getMetaInfo();
         // initial spacing and position
         float spacing[3]{.0,.0,.0};
-        spacing[0] = getSpacing(startPoint,endPoint,imageXYZ,0);
+        spacing[0] = getSpacing(startPoint,endPoint,aOriginXYZ,0);
-        spacing[1] = getSpacing(startPoint,endPoint,imageXYZ,1);
+        spacing[1] = getSpacing(startPoint,endPoint,aOriginXYZ,1);
-        spacing[2] = getSpacing(startPoint,endPoint,imageXYZ,2);
+        spacing[2] = getSpacing(startPoint,endPoint,aOriginXYZ,2);
-
+        float originPosition[3];
-        float originPosition[3] = {
+        initOriginPosition(type, originPosition, endPoint, aOriginXYZ, spacing);
-            endPoint[1]*(float)1000.0,
+        initPixelSpacing(type, spacing, Rows, Cols);
-            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++)
@@ -234,10 +390,6 @@ namespace Recon
         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;
@@ -246,19 +398,29 @@ namespace Recon
             windowCenter = (windowCenter-min)/dSlope;
             windowWidth = (windowWidth-min)/dSlope;
         }
+        initImageOrientation(type);
+        dcmGenerateUniqueIdentifier(mSeriesInstanceUID,SITE_SERIES_UID_ROOT);
+        dataset->putAndInsertString(DCM_SeriesInstanceUID, mSeriesInstanceUID);
         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;
+        ushort* udata = new ushort[sliceSize]{0};
         for (size_t i = 0; i < Slices; i++)
         {
+            Aurora::Matrix image = Aurora::Matrix::copyFromRawData(data + i * sliceSize, Rows, Cols);
+            Aurora::Matrix dicomImage = transferToDicomType(image, type);
+            for (size_t i = 0; i < sliceSize; i++)
+            {
+                udata[i] = (ushort)((dicomImage[i] - min)/dSlope);
+            }
             dataset->putAndInsertString(DCM_AccessionNumber, mPatientData.getAccessionNumber().data());
             dataset->putAndInsertString(DCM_InstanceNumber, to_string(i).data());
-            string pos = to_string(originPosition[1] ) + "\\" +
+            string pos = to_string(originPosition[0] ) + "\\" +
-                        to_string(originPosition[2]+ i * spacing[2]) + "\\" +
+            to_string(originPosition[1]+ i * spacing[2]) + "\\" +
-                        to_string(originPosition[0]);
+            to_string(originPosition[2]);
             dataset->putAndInsertString(DCM_ImagePositionPatient, pos.data());
             string loc = to_string(originLocation + i * spacing[2]);
             dataset->putAndInsertString(DCM_SliceLocation, loc.data());
@@ -269,12 +431,12 @@ namespace Recon
             metaInfo->putAndInsertString(DCM_MediaStorageSOPInstanceUID, SOP_UID);
             dataset->putAndInsertString(DCM_SOPInstanceUID, SOP_UID);
             
-            dataset->putAndInsertUint16Array(DCM_PixelData, udata + sliceSize * i,
+            dataset->putAndInsertUint16Array(DCM_PixelData, udata,
                                             sliceSize);
             string path2 = path+to_string(i)+".dcm";
             OFFilename file ;
             file = path2.data();
-            auto status =dcmFile.saveFile(file, EXS_LittleEndianExplicit);
+            auto status =mDicomFile->saveFile(file, EXS_LittleEndianExplicit);
             if(status.bad())
             {
                 RECON_ERROR("Save dicom file {0} fail, beacuse of {1}\r\n",path2.data(),status.text());
@@ -283,5 +445,6 @@ namespace Recon
                 RECON_INFO("Save dicom file {0} \r\n",path2.data());
             }
         }
+        delete udata;
     }
 }

src/common/DICOMExporter.h

@@ -7,6 +7,10 @@
 #include "Parser/Data/PatientData.h"
 #include "Parser/Data/MetaData.h"
 
+class DcmDataset;
+class DcmMetaInfo;
+class DcmFileFormat;
+
 namespace Recon {
 
 class DICOMExporter
@@ -15,7 +19,10 @@ public:
     enum ImageType{
         REFL,
         SOS,
-        ATT
+        ATT,
+        CoronalPositioning,
+        SagittalPositioning,
+        TransversePositioning
     };
     explicit DICOMExporter(const PatientData& aPatientData, const MetaData& aMetaData);
     //仅仅用于测试！！！
@@ -27,14 +34,25 @@ public:
     ~DICOMExporter();
     void setExportBasePath(const std::string & path);
     void exportDICOM(Aurora::Matrix aMatrix, ImageType type);
+    void exportDICOM(Aurora::Matrix aMatrix, const Aurora::Matrix& aOriginMatrix, ImageType type);
+
 private:
     std::string mBasePath;
     std::string mStudyDate;
     std::string mStudyTime;
     std::string mStudyInstanceUID;
+    std::string mLaterality;
+    char mFrameOfReferenceUID[100]={0};
     char mSeriesInstanceUID[100]={0};
     PatientData mPatientData;
-    void exportDCM(std::string path, float * startPoint, float* endPoint, int* imageXYZ,float* data, int type);
+    DcmFileFormat* mDicomFile;
+    void exportDCM(std::string path, float * startPoint, float* endPoint, int* imageXYZ,float* data, ImageType type, int* aOriginXYZ);
+    std::string getImageOrientationPatient(DICOMExporter::ImageType aType);
+    void initialDicomFile();
+    void initImageOrientation(ImageType aType);
+    void initOriginPosition(ImageType aType, float* aOriginPosition, float* aPoint, int* aXYZ, float* aSpacing);
+    void initPixelSpacing(ImageType aType, float* aSpacing, int aRow, int aColumn);
+
 };