diff --git a/.gitignore b/.gitignore index 24ccbf0b..39d1724b 100644 --- a/.gitignore +++ b/.gitignore @@ -35,3 +35,5 @@ CMakeLists.txt.user* __pycache__/ *.py[cod] *$py.class + +node_modules diff --git a/README.md b/README.md index 868b0843..394703d9 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,4 @@ - -[![OpenHub](https://www.openhub.net/p/dcmqi/widgets/project_thin_badge.gif)](https://www.openhub.net/p/dcmqi) [![codecov](https://codecov.io/gh/QIICR/dcmqi/branch/master/graph/badge.svg)](https://codecov.io/gh/QIICR/dcmqi) [![Join the chat at https://gitter.im/QIICR/dcmqi](https://badges.gitter.im/QIICR/dcmqi.svg)](https://gitter.im/QIICR/dcmqi?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge) - - +[![OpenHub](https://www.openhub.net/p/dcmqi/widgets/project_thin_badge.gif)](https://www.openhub.net/p/dcmqi) [![codecov](https://codecov.io/gh/QIICR/dcmqi/branch/master/graph/badge.svg)](https://codecov.io/gh/QIICR/dcmqi) [![Join the chat at https://gitter.im/QIICR/dcmqi](https://badges.gitter.im/QIICR/dcmqi.svg)](https://gitter.im/QIICR/dcmqi?utm_source=badge&utm_medium=badge&utm_campaign=pr-badge&utm_content=badge) | Docker | [![](https://images.microbadger.com/badges/version/qiicr/dcmqi:v1.0.5.svg)](https://microbadger.com/images/qiicr/dcmqi:v1.0.5) | [![](https://images.microbadger.com/badges/version/qiicr/dcmqi.svg)](https://microbadger.com/images/qiicr/dcmqi) | |--------|--------------------------------------------------------------------------------------------------------------------------------|-------------------------------------------------------------------------------------------------------------------| diff --git a/apps/paramaps/itkimage2paramap.cxx b/apps/paramaps/itkimage2paramap.cxx index 583b86d9..22f5b6f9 100644 --- a/apps/paramaps/itkimage2paramap.cxx +++ b/apps/paramaps/itkimage2paramap.cxx @@ -3,7 +3,7 @@ // DCMQI includes #undef HAVE_SSTREAM // Avoid redefinition warning -#include "dcmqi/ParaMapConverter.h" +#include "dcmqi/ParametricMapConverter.h" #include "dcmqi/internal/VersionConfigure.h" @@ -44,7 +44,7 @@ int main(int argc, char *argv[]) std::string metadata( (std::istreambuf_iterator(metainfoStream) ), (std::istreambuf_iterator())); - DcmDataset* result = dcmqi::ParaMapConverter::itkimage2paramap(parametricMapImage, dcmDatasets, metadata); + DcmDataset* result = dcmqi::itkimage2paramapReplacement(parametricMapImage, dcmDatasets, metadata); if (result == NULL) { return EXIT_FAILURE; diff --git a/apps/paramaps/paramap2itkimage.cxx b/apps/paramaps/paramap2itkimage.cxx index aa1def1d..13161cbb 100644 --- a/apps/paramaps/paramap2itkimage.cxx +++ b/apps/paramaps/paramap2itkimage.cxx @@ -3,7 +3,7 @@ // DCMQI includes #undef HAVE_SSTREAM // Avoid redefinition warning -#include "dcmqi/ParaMapConverter.h" +#include "dcmqi/ParametricMapConverter.h" #include "dcmqi/internal/VersionConfigure.h" @@ -25,7 +25,7 @@ int main(int argc, char *argv[]) CHECK_COND(sliceFF.loadFile(inputFileName.c_str())); DcmDataset* dataset = sliceFF.getDataset(); - pair result = dcmqi::ParaMapConverter::paramap2itkimage(dataset); + pair result = dcmqi::paramap2itkimageReplacement(dataset); string fileExtension = helper::getFileExtensionFromType(outputType); diff --git a/apps/seg/itkimage2segimage.cxx b/apps/seg/itkimage2segimage.cxx index 4619554e..cd2c912e 100644 --- a/apps/seg/itkimage2segimage.cxx +++ b/apps/seg/itkimage2segimage.cxx @@ -3,7 +3,7 @@ // DCMQI includes #undef HAVE_SSTREAM // Avoid redefinition warning -#include "dcmqi/ImageSEGConverter.h" +#include "dcmqi/SegmentationImageConverter.h" #include "dcmqi/internal/VersionConfigure.h" typedef dcmqi::Helper helper; @@ -91,7 +91,7 @@ int main(int argc, char *argv[]) segmentations = segmentationsReordered; } - DcmDataset* result = dcmqi::ImageSEGConverter::itkimage2dcmSegmentation(dcmDatasets, segmentations, metadata, skipEmptySlices); + DcmDataset* result = dcmqi::SegmentationImageConverter::itkimage2dcmSegmentation(dcmDatasets, segmentations, metadata, skipEmptySlices); if (result == NULL){ return EXIT_FAILURE; diff --git a/apps/seg/segimage2itkimage.cxx b/apps/seg/segimage2itkimage.cxx index 7784e311..82721a36 100644 --- a/apps/seg/segimage2itkimage.cxx +++ b/apps/seg/segimage2itkimage.cxx @@ -3,7 +3,7 @@ // DCMQI includes #undef HAVE_SSTREAM // Avoid redefinition warning -#include "dcmqi/ImageSEGConverter.h" +#include "dcmqi/SegmentationImageConverter.h" #include "dcmqi/internal/VersionConfigure.h" @@ -24,7 +24,7 @@ int main(int argc, char *argv[]) CHECK_COND(sliceFF.loadFile(inputSEGFileName.c_str())); DcmDataset* dataset = sliceFF.getDataset(); - pair , string> result = dcmqi::ImageSEGConverter::dcmSegmentation2itkimage(dataset); + pair , string> result = dcmqi::dcmSegmentation2itkimageReplacement(dataset); string outputPrefix = prefix.empty() ? "" : prefix + "-"; diff --git a/doc/examples/pm-example-float.json b/doc/examples/pm-example-float.json index 7e70af1b..ab3f3528 100644 --- a/doc/examples/pm-example-float.json +++ b/doc/examples/pm-example-float.json @@ -22,8 +22,8 @@ "CodeMeaning": "m2/s" }, "MeasurementMethodCode": { - "CodeValue": "DWMPxxxx10", - "CodingSchemeDesignator": "99QIICR", + "CodeValue": "113290", + "CodingSchemeDesignator": "DCM", "CodeMeaning": "Mono-exponential diffusion model" }, "SourceImageDiffusionBValues": ["0","1400"], diff --git a/doc/examples/pm-example.json b/doc/examples/pm-example.json index f0b2970d..30df68f1 100644 --- a/doc/examples/pm-example.json +++ b/doc/examples/pm-example.json @@ -22,8 +22,8 @@ "CodeMeaning": "um2/s" }, "MeasurementMethodCode": { - "CodeValue": "DWMPxxxx10", - "CodingSchemeDesignator": "99QIICR", + "CodeValue": "113290", + "CodingSchemeDesignator": "DCM", "CodeMeaning": "Mono-exponential diffusion model" }, "SourceImageDiffusionBValues": ["0","1400"], diff --git a/include/dcmqi/DICOMFrame.h b/include/dcmqi/DICOMFrame.h new file mode 100644 index 00000000..a2ef39d3 --- /dev/null +++ b/include/dcmqi/DICOMFrame.h @@ -0,0 +1,104 @@ +// +// Created by Andrey Fedorov on 3/9/17. +// + +#ifndef DCMQI_DICOMFRAME_H +#define DCMQI_DICOMFRAME_H + +#include +#include +#include "ImageVolumeGeometry.h" +#include "Exceptions.h" + +namespace dcmqi { + + // Describe input/output frames for the purposes of sorting and associating with the + // slices of the ITK volume + class DICOMFrame { + public: + // distinguish between the frames from the legacy data and enhanced multiframe objects + enum { + LegacyInstanceFrame = 0, + EnhancedInstanceFrame + }; + + DICOMFrame(DcmDataset *dataset, int number=0) : + frameNumber(number), + frameDataset(dataset) { + + Uint32 numFrames; + if(dataset->findAndGetUint32(DCM_NumberOfFrames, numFrames).good()){ + frameType = EnhancedInstanceFrame; + if(!number){ + std::cerr << "ERROR: DICOMFrame object for an enhanced frame is initialized with frame 0!" << std::endl; + } + } else { + frameType = LegacyInstanceFrame; + initializeFrameGeometryFromLegacyInstance(); + } + + OFString seriesUIDOF, instanceUIDOF, classUIDOF; + if(dataset->findAndGetOFString(DCM_SeriesInstanceUID, seriesUIDOF).good()){ + seriesUID = seriesUIDOF.c_str(); + } + if(dataset->findAndGetOFString(DCM_SOPInstanceUID, instanceUIDOF).good()){ + instanceUID = instanceUIDOF.c_str(); + } + if(dataset->findAndGetOFString(DCM_SOPClassUID, classUIDOF).good()){ + classUID = classUIDOF.c_str(); + } + + }; + + int getFrameNumber() const { + return frameNumber; // 0 for legacy datasets, 1 or above for enhanced objects + }; + + vnl_vector getFrameIPP(){ + return frameIPP; + }; + + string getSeriesUID() const { + return seriesUID; + } + + string getInstanceUID() const{ + return instanceUID; + } + + string getClassUID() const { + return classUID; + } + + DcmDataset* getDataset() const { + return frameDataset; + } + + private: + + int initializeFrameGeometryFromLegacyInstance(); + int initializeFrameGeometryFromEnhancedInstance(); + + int frameType; + DcmDataset *frameDataset; + int frameNumber; + vnl_vector frameIPP; + + string seriesUID, instanceUID, classUID; + + ImageVolumeGeometry frameGeometry; + + }; + + struct DICOMFrame_compare { + bool operator() (const DICOMFrame& lhs, const DICOMFrame& rhs) const{ + std::stringstream s1,s2; + s1 << lhs.getInstanceUID(); + s2 << rhs.getInstanceUID(); + return (s1.str() < s2.str()) && (lhs.getFrameNumber() < rhs.getFrameNumber()); + } + }; + +} + +#endif //DCMQI_DICOMFRAME_H diff --git a/include/dcmqi/Helper.h b/include/dcmqi/Helper.h index 903bd8d0..483d773d 100644 --- a/include/dcmqi/Helper.h +++ b/include/dcmqi/Helper.h @@ -17,8 +17,9 @@ // DCMQI includes #include "dcmqi/Exceptions.h" -using namespace std; +#include +using namespace std; namespace dcmqi { @@ -60,10 +61,21 @@ namespace dcmqi { static CodeSequenceMacro stringToCodeSequenceMacro(string str); static DSRCodedEntryValue stringToDSRCodedEntryValue(string str); + static string codeSequenceMacroToString(CodeSequenceMacro); + + static CodeSequenceMacro jsonToCodeSequenceMacro(Json::Value); static void checkValidityOfFirstSrcImage(DcmSegmentation *segdoc); static CodeSequenceMacro* createNewCodeSequence(const string& code, const string& designator, const string& meaning); + + static OFString generateUID(); + static OFString getTagAsOFString(DcmDataset*, DcmTagKey); + + static string getCodeSequenceValue(CodeSequenceMacro &codeSequence); + static string getCodeSequenceDesignator(CodeSequenceMacro &codeSequence); + static string getCodeSequenceMeaning(CodeSequenceMacro &codeSequence); + static Json::Value codeSequence2Json(CodeSequenceMacro &codeSequence); }; } diff --git a/include/dcmqi/ImageVolume.h b/include/dcmqi/ImageVolume.h new file mode 100644 index 00000000..b6ebf83e --- /dev/null +++ b/include/dcmqi/ImageVolume.h @@ -0,0 +1,74 @@ +// +// Created by Andrey Fedorov on 3/9/17. +// + +#ifndef DCMQI_IMAGEVOLUME_H +#define DCMQI_IMAGEVOLUME_H + +#include +#include +#include + +#include + +namespace dcmqi { + + // Maintain properties of the image volume + // Attributes ara parallel to those of itkImageData, but limited to what we need for the task of conversion, + // and the class is not templated over the pixel type, since we may not know the pixel type + // at the time class is instantiated. + // + // Initially, limit implementation and support to just Float32 used by the original PM converter. + + class ImageVolume { + public: + // pixel types that are relevant for the types of objects we want to support + enum { + FLOAT32 = 0, // PM + FLOAT64, // PM + UINT16 // PM or SEG + }; + + typedef IODFloatingPointImagePixelModule::value_type Float32PixelType; + typedef itk::Image Float32ITKImageType; + + ImageVolume(){ + rowDirection.set_size(3); + columnDirection.set_size(3); + sliceDirection.set_size(3); + origin.set_size(3); + spacing.set_size(3); + pixelData = NULL; + } + + // while going from DICOM PM/SEG, we get volume information from FGInterface + int initializeFromDICOM(FGInterface&); + int initializeFromITK(Float32ITKImageType::Pointer); + + // TODO - inherited class? or separate segments before passing to this one? + // int initializeFromSegment(FGInterface&, unsigned); + + protected: + int initializeDirections(FGInterface &); + int initializeExtent(FGInterface &); + bool getDeclaredSliceSpacing(FGInterface&); + bool getCalculatedSliceSpacing(); + + int setDirections(vnl_vector rowDirection, vnl_vector columnDirection, vnl_vector sliceDirection); + int setOrigin(vnl_vector); + + private: + + // use vnl_vector to simplify support of vector calculations + vnl_vector rowDirection, columnDirection, sliceDirection; + vnl_vector origin; + unsigned sliceExtent; + vnl_vector spacing; + void* pixelData; + int pixelDataType; + }; + +}; + + +#endif //DCMQI_IMAGEVOLUME_H diff --git a/include/dcmqi/ImageVolumeGeometry.h b/include/dcmqi/ImageVolumeGeometry.h new file mode 100644 index 00000000..c437cc89 --- /dev/null +++ b/include/dcmqi/ImageVolumeGeometry.h @@ -0,0 +1,98 @@ +// +// Created by Andrey Fedorov on 3/11/17. +// + +#ifndef DCMQI_IMAGEVOLUMEGEOMETRY_H +#define DCMQI_IMAGEVOLUMEGEOMETRY_H + + +#include +#include +#include +#include +#include + +class ImageVolumeGeometry { + + friend class MultiframeObject; + friend class ParametricMapObject; + friend class SegmentationImageObject; + +public: + typedef itk::Vector DoubleVectorType; + typedef itk::Size<3> SizeType; + + // implementation of the image volume geometry + // NB: duplicated from MultiframeObject! + typedef unsigned char DummyPixelType; + typedef itk::Image DummyImageType; + typedef DummyImageType::PointType PointType; + typedef DummyImageType::DirectionType DirectionType; + + ImageVolumeGeometry(); + // initialize from DICOM + ImageVolumeGeometry(DcmDataset*); + + int setSpacing(DoubleVectorType); + int setOrigin(PointType); + int setExtent(SizeType); + int setDirections(DirectionType); + + template + typename T::Pointer getITKRepresentation(){ + typename T::Pointer image; + typename T::IndexType index; + typename T::SizeType size; + typename T::DirectionType direction; + typename T::SpacingType spacing; + typename T::RegionType region; + typename T::PointType origin; + + image = T::New(); + + index.Fill(0); + + size[0] = extent[0]; + size[1] = extent[1]; + size[2] = extent[2]; + + region.SetIndex(index); + region.SetSize(size); + + spacing[0] = this->spacing[0]; + spacing[1] = this->spacing[1]; + spacing[2] = this->spacing[2]; + + for (int i = 0; i < 3; i++) + direction[i][0] = rowDirection[i]; + for (int i = 0; i < 3; i++) + direction[i][1] = columnDirection[i]; + for (int i = 0; i < 3; i++) + direction[i][2] = sliceDirection[i]; + + origin[0] = this->origin[0]; + origin[1] = this->origin[1]; + origin[2] = this->origin[2]; + + image->SetRegions(region); + image->SetOrigin(origin); + image->SetDirection(direction); + image->SetSpacing(spacing); + + image->Allocate(); + image->FillBuffer(0); + + return image; + } + +protected: + // use vnl_vector to simplify support of vector calculations + vnl_vector rowDirection, columnDirection, sliceDirection; + vnl_vector origin; + vnl_vector extent; + vnl_vector spacing; + +}; + + +#endif //DCMQI_IMAGEVOLUMEGEOMETRY_H diff --git a/include/dcmqi/JSONSegmentationMetaInformationHandler.h b/include/dcmqi/JSONSegmentationMetaInformationHandler.h index 219b0087..7a893b15 100644 --- a/include/dcmqi/JSONSegmentationMetaInformationHandler.h +++ b/include/dcmqi/JSONSegmentationMetaInformationHandler.h @@ -47,6 +47,7 @@ namespace dcmqi { void readSegmentAttributes(); Json::Value createAndGetSegmentAttributes(); + }; } diff --git a/include/dcmqi/ConverterBase.h b/include/dcmqi/MultiframeConverter.h similarity index 95% rename from include/dcmqi/ConverterBase.h rename to include/dcmqi/MultiframeConverter.h index ad7e9f49..84089ab9 100644 --- a/include/dcmqi/ConverterBase.h +++ b/include/dcmqi/MultiframeConverter.h @@ -35,15 +35,21 @@ using namespace std; +// common type definitions typedef short ShortPixelType; typedef itk::Image ShortImageType; typedef itk::ImageFileReader ShortReaderType; namespace dcmqi { - class ConverterBase { + class MultiframeConverter { + public: + //virtual int convert(); + MultiframeConverter(){ + }; + + - protected: static IODGeneralEquipmentModule::EquipmentInfo getEquipmentInfo(); static IODEnhGeneralEquipmentModule::EquipmentInfo getEnhEquipmentInfo(); static ContentIdentificationMacro createContentIdentificationInformation(JSONMetaInformationHandlerBase &metaInfo); @@ -124,7 +130,7 @@ namespace dcmqi { */ // Determine ordering of the frames, keep mapping from ImagePositionPatient string - // to the distance, and keep track (just out of curiousity) how many frames overlap + // to the distance, and keep track (just out of curiosity) how many frames overlap vnl_vector refOrigin(3); { OFBool isPerFrame; @@ -279,6 +285,10 @@ namespace dcmqi { CHECK_COND(dcmDatasets[i]->findAndGetOFString(DCM_ImagePositionPatient, ippStr, j)); ippPoint[j] = atof(ippStr.c_str()); } + // NB: this will map slice origin to index without failure, unless the point is out + // of FOV bounds! + // TODO: do a better job matching volume slices by considering comparison of the origin + // and orientation of the slice within tolerance if(!labelImage->TransformPhysicalPointToIndex(ippPoint, ippIndex)){ //cout << "image position: " << ippPoint << endl; //cerr << "ippIndex: " << ippIndex << endl; diff --git a/include/dcmqi/MultiframeObject.h b/include/dcmqi/MultiframeObject.h new file mode 100644 index 00000000..96d85632 --- /dev/null +++ b/include/dcmqi/MultiframeObject.h @@ -0,0 +1,185 @@ +// +// Created by Andrey Fedorov on 3/11/17. +// + +#ifndef DCMQI_MULTIFRAMEOBJECT_H +#define DCMQI_MULTIFRAMEOBJECT_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include "vnl/vnl_cross.h" + +#include +#include +#include +#include +#include + +#include "dcmqi/Exceptions.h" +#include "dcmqi/ImageVolumeGeometry.h" +#include "dcmqi/DICOMFrame.h" + +using namespace std; +/* + * Class to support conversion between DICOM and ITK representations + * of multiframe objects + * + * + */ +class MultiframeObject { + +public: + + // initialization from DICOM is always from the dataset(s) + // that encode DICOM representation, optionally augmented by the + // derivation datasets that potentially can be helpful in some + // situations. + int initializeFromDICOM(std::vector sourceDataset); + int initializeFromDICOM(std::vector sourceDatasets, + std::vector derivationDatasets); + + // initialization from ITK will be specific to the derived types, + // although there will be some common initialization of the metadata + + // Output is always a single DcmDataset, since this is a multiframe + // object + DcmDataset* getDcmDataset() const; + Json::Value getMetaDataJson() const { + return metaDataJson; + }; + + // get ITK representation will be specific to the derived classes, + // since the type of the ITK image and the number of ITK images is + // different between PM and SEG + +protected: + + // Helpers to convert to dummy image type to support common + // implementation of the image volume geometry + typedef unsigned char DummyPixelType; + typedef itk::Image DummyImageType; + typedef DummyImageType::PointType PointType; + typedef DummyImageType::DirectionType DirectionType; + typedef DummyImageType::SpacingType SpacingType; + typedef DummyImageType::SizeType SizeType; + + int initializeMetaDataFromString(const std::string&); + // what this function does depends on whether we are coming from + // DICOM or from ITK. No parameters, since all it does is exchange + // between DICOM and MetaData + int initializeContentIdentification(); + + // from ITK + int initializeVolumeGeometryFromITK(DummyImageType::Pointer); + + int initializeVolumeGeometryFromDICOM(FGInterface &fgInterface); + int getImageDirections(FGInterface& fgInterface, DirectionType &dir); + + int computeVolumeExtent(FGInterface& fgInterface, vnl_vector &sliceDirection, PointType &imageOrigin, + double &sliceSpacing, double &volumeExtent); + int getDeclaredImageSpacing(FGInterface &fgInterface, SpacingType &spacing); + + // initialize attributes of the composite context that are common for all multiframe objects + int initializeCompositeContext(); + int initializeSeriesSpecificAttributes(IODGeneralSeriesModule&, IODGeneralImageModule&); + // check whether all of the attributes required for initialization of the object are present in the + // input metadata + virtual bool metaDataIsComplete(); + + ContentItemMacro* initializeContentItemMacro(CodeSequenceMacro, CodeSequenceMacro); + + // List of tags, and FGs they belong to, for initializing dimensions module + int initializeDimensions(std::vector >); + int initializePixelMeasuresFG(); + int initializePlaneOrientationFG(); + int initializeCommonInstanceReferenceModule(IODCommonInstanceReferenceModule &, vector >&); + + int mapVolumeSlicesToDICOMFrames(const vector, + vector >&); + + static std::vector findIntersectingSlices(ImageVolumeGeometry& volume, dcmqi::DICOMFrame& frame); + + int addDerivationItemToDerivationFG(FGDerivationImage* fgder, set frames, + CodeSequenceMacro purposeOfReferenceCode = CodeSequenceMacro("121322","DCM","Source image for image processing operation"), + CodeSequenceMacro derivationCode = CodeSequenceMacro("110001","DCM","Image Processing")); + + void insertDerivationSeriesInstance(string seriesUID, string instanceUID); + + int findIntersectingSlicesAndAddDerivationSeriesInstance( + vector > &slice2frame, DcmDataset *dcm, int frameNo=0); + + int setDerivationDatasets(std::vector derivationDatasets){ + for(std::vector::const_iterator vIt=derivationDatasets.begin(); + vIt!=derivationDatasets.end();++vIt) + derivationDcmDatasets.push_back(*vIt); + return EXIT_SUCCESS; + } + + // constants to describe original representation of the data being converted + enum { + DICOM_REPR = 0, + ITK_REPR + }; + int sourceRepresentationType; + + // TODO: abstract this into a different class, which would help with: + // - checking for presence of attributes + // - handling of defaults (in the future, initialized from the schema?) + // - simplifying common access patterns (access to the Code tuples) + Json::Value metaDataJson; + + // Multiframe DICOM object representation + // probably not needed, since need object-specific DCMTK class in + // derived classes + DcmDataset* dcmRepresentation; + + // probably not needed at this level, since for SEG each segment will + // have separate geometry definition + ImageVolumeGeometry volumeGeometry; + + // DcmDataset(s) that hold the original representation of the + // object, when the sourceRepresentationType == DICOM_REPR + OFVector sourceDcmDatasets; + + // Common components present in the derived classes + ContentIdentificationMacro contentIdentificationMacro; + IODMultiframeDimensionModule dimensionsModule; + + // DcmDataset(s) that were used to derive the object + // Probably will only be populated when sourceRepresentationType == ITK_REPR + // Purpose of those: + // 1) initialize derivation derivationImageFG (ITK->) + // 2) initialize CommonInstanceReferenceModule (ITK->) + // 3) initialize common attributes (ITK->) + OFVector derivationDcmDatasets; + + // Functional groups common to all MF objects: + // - Shared + FGPixelMeasures pixelMeasuresFG; + FGPlaneOrientationPatient planeOrientationPatientFG; + // - Per-frame + OFVector planePosPatientFGList; + OFVector frameContentFGList; + OFVector derivationImageFGList; + + // Mapping from the derivation items SeriesUIDs to InstanceUIDs + std::map > derivationSeriesToInstanceUIDs; + + vnl_vector getFrameOrigin(FGInterface &fgInterface, int frameId) const; + vnl_vector getFrameOrigin(FGPlanePosPatient *planposfg) const; +}; + + +#endif //DCMQI_MULTIFRAMEOBJECT_H diff --git a/include/dcmqi/ParaMapConverter.h b/include/dcmqi/ParametricMapConverter.h similarity index 52% rename from include/dcmqi/ParaMapConverter.h rename to include/dcmqi/ParametricMapConverter.h index fcabbf49..57c244d3 100644 --- a/include/dcmqi/ParaMapConverter.h +++ b/include/dcmqi/ParametricMapConverter.h @@ -18,8 +18,9 @@ #include // DCMQI includes -#include "dcmqi/ConverterBase.h" +#include "dcmqi/MultiframeConverter.h" #include "dcmqi/JSONParametricMapMetaInformationHandler.h" +#include "dcmqi/ImageVolume.h" typedef IODFloatingPointImagePixelModule::value_type FloatPixelType; typedef itk::Image FloatImageType; @@ -28,22 +29,53 @@ typedef itk::MinimumMaximumImageCalculator MinMaxCalculatorType; using namespace std; - namespace dcmqi { - class ParaMapConverter : public ConverterBase { + DcmDataset* itkimage2paramapReplacement(const FloatImageType::Pointer ¶metricMapImage, vector dcmDatasets, + const string &metaData); + pair paramap2itkimageReplacement(DcmDataset *pmapDataset); + + class ParametricMapConverter : public MultiframeConverter { public: - static DcmDataset* itkimage2paramap(const FloatImageType::Pointer ¶metricMapImage, vector dcmDatasets, - const string &metaData); + // placeholder to replace static function going from ITK to DICOM + ParametricMapConverter(const FloatImageType::Pointer ¶metricMapImage, vector dcmDatasets, + const string &metaData); + + // placeholder to replace static function going from DICOM to ITK + ParametricMapConverter(DcmDataset*); + + static DcmDataset* itkimage2paramap(const FloatImageType::Pointer ¶metricMapImage, vector dcmDatasets, + const string &metaData); static pair paramap2itkimage(DcmDataset *pmapDataset); + + // given one representation, generate the parallel one + int convert(); + + // get the result + FloatImageType::Pointer getFloat32ITKImage() const; + string getMetaData() const; + DcmDataset* getDataset() const; + protected: + + ParametricMapConverter() : parametricMapVolume(NULL), parametricMapDataset(NULL) {}; + static OFCondition addFrame(DPMParametricMapIOD &map, const FloatImageType::Pointer ¶metricMapImage, const JSONParametricMapMetaInformationHandler &metaInfo, const unsigned long frameNo, OFVector perFrameGroups); static void populateMetaInformationFromDICOM(DcmDataset *pmapDataset, JSONParametricMapMetaInformationHandler &metaInfo); + private: + // these are the items we convert to and from, essentially + ImageVolume* parametricMapVolume; + DcmDataset* parametricMapDataset; + + // these are the items we will need in the process of conversion + vector referencedDatasets; + string metaData; + }; } diff --git a/include/dcmqi/ParametricMapObject.h b/include/dcmqi/ParametricMapObject.h new file mode 100644 index 00000000..b34e34ea --- /dev/null +++ b/include/dcmqi/ParametricMapObject.h @@ -0,0 +1,92 @@ +// +// Created by Andrey Fedorov on 3/11/17. +// + +#ifndef DCMQI_PARAMETRICMAPOBJECT_H +#define DCMQI_PARAMETRICMAPOBJECT_H + + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "MultiframeObject.h" + +/* + * + */ +class ParametricMapObject : public MultiframeObject { +public: + + ParametricMapObject(){ + parametricMap = NULL; + } + + typedef IODFloatingPointImagePixelModule::value_type Float32PixelType; + typedef itk::Image Float32ITKImageType; + + // metadata is mandatory, since not all of the attributes can be present + // in the derivation DcmDataset(s) + int initializeFromITK(Float32ITKImageType::Pointer, const string&); + // metadata is mandatory, optionally, derivation DcmDataset(s) can + // help + int initializeFromITK(Float32ITKImageType::Pointer, const string&, + std::vector); + + int updateMetaDataFromDICOM(std::vector); + + int getDICOMRepresentation(DcmDataset& dcm){ + if(parametricMap) + CHECK_COND(parametricMap->write(dcm)); + }; + + int initializeFromDICOM(DcmDataset * sourceDataset); + + void initializeMetaDataFromDICOM(); + + Float32ITKImageType::Pointer getITKRepresentation() const { + return itkImage; + } + +protected: + typedef itk::CastImageFilter + Float32ToDummyCasterType; + typedef itk::MinimumMaximumImageCalculator MinMaxCalculatorType; + + int initializeEquipmentInfo(); + int initializeVolumeGeometry(); + int createDICOMParametricMap(); + int createITKParametricMap(); + int createITKImageFromFrames(FGInterface&, DPMParametricMapIOD::Frames); + int initializeCompositeContext(); + int initializeFrameAnatomyFG(); + int initializeRWVMFG(); + int initializeFrames(vector >&); + + bool isDerivationFGRequired(vector >& slice2frame); + + IODEnhGeneralEquipmentModule::EquipmentInfo enhancedEquipmentInfoModule; + + // Functional groups specific to PM: + // - Shared + FGFrameAnatomy frameAnatomyFG; + FGIdentityPixelValueTransformation identityPixelValueTransformationFG; + FGParametricMapFrameType parametricMapFrameTypeFG; + FGRealWorldValueMapping rwvmFG; + + // Data containers specific to this object + Float32ITKImageType::Pointer itkImage; + +private: + DPMParametricMapIOD* parametricMap; + CodeSequenceMacro *createCodeSequenceFromMetadata(const string &codeName) const; +}; + + +#endif //DCMQI_PARAMETRICMAPOBJECT_H diff --git a/include/dcmqi/SegmentVolume.h b/include/dcmqi/SegmentVolume.h new file mode 100644 index 00000000..504693cd --- /dev/null +++ b/include/dcmqi/SegmentVolume.h @@ -0,0 +1,17 @@ +// +// Created by Andrey Fedorov on 3/9/17. +// + +#ifndef DCMQI_SEGMENTVOLUME_H +#define DCMQI_SEGMENTVOLUME_H + +#include "dcmqi/ImageVolume.h" + +namespace dcmqi { + class SegmentVolume : public ImageVolume { + + }; +} + + +#endif //DCMQI_SEGMENTVOLUME_H diff --git a/include/dcmqi/ImageSEGConverter.h b/include/dcmqi/SegmentationImageConverter.h similarity index 77% rename from include/dcmqi/ImageSEGConverter.h rename to include/dcmqi/SegmentationImageConverter.h index deb8e7f4..75c9d506 100644 --- a/include/dcmqi/ImageSEGConverter.h +++ b/include/dcmqi/SegmentationImageConverter.h @@ -22,7 +22,7 @@ #include // DCMQI includes -#include "dcmqi/ConverterBase.h" +#include "dcmqi/MultiframeConverter.h" #include "dcmqi/JSONSegmentationMetaInformationHandler.h" using namespace std; @@ -32,7 +32,12 @@ typedef itk::LabelImageToLabelMapFilter LabelToLabelMapFilterTyp namespace dcmqi { - class ImageSEGConverter : public ConverterBase { + DcmDataset* itkimage2paramapReplacement(vector dcmDatasets, vector segmentations, + const string &metaData, bool skipEmptySlices); + + pair , string> dcmSegmentation2itkimageReplacement(DcmDataset *segDataset); + + class SegmentationImageConverter : public MultiframeConverter { public: static DcmDataset* itkimage2dcmSegmentation(vector dcmDatasets, diff --git a/include/dcmqi/SegmentationImageObject.h b/include/dcmqi/SegmentationImageObject.h new file mode 100644 index 00000000..53a699aa --- /dev/null +++ b/include/dcmqi/SegmentationImageObject.h @@ -0,0 +1,88 @@ +// +// Created by Andrey Fedorov on 3/11/17. +// + +#ifndef DCMQI_SEGMENTATIONIMAGEOBJECT_H +#define DCMQI_SEGMENTATIONIMAGEOBJECT_H + +// DCMQI includes +#include "MultiframeObject.h" +#include "Helper.h" +#include "QIICRConstants.h" + +// DCMTK includes +#include +#include +#include +#include +#include +#include + +// ITK includes +#include +#include + + +class SegmentationImageObject : public MultiframeObject { + +public: + + typedef short ShortPixelType; + typedef itk::Image ShortImageType; + + SegmentationImageObject(){ + segmentation = NULL; + } + + int initializeFromITK(vector, const string &, vector, bool); + + int initializeFromDICOM(DcmDataset* sourceDataset); + + int getDICOMRepresentation(DcmDataset& dcm){ + if(segmentation) + CHECK_COND(segmentation->writeDataset(dcm)); + }; + map getITKRepresentation() const { + // TODO: think about naming + return segment2image; + } + +protected: + + typedef itk::CastImageFilter ShortToDummyCasterType; + // Data containers specific to this object + ShortImageType::Pointer itkImage; + vector itkImages; + + // ITK images corresponding to the individual segments + map segment2image; + + DcmSegmentation* segmentation; + + int initializeEquipmentInfo(); + int initializeVolumeGeometry(); + int initializeCompositeContext(); + + // returns a vector with a size equal to the number of frames each holding segmentID and sliceNumber + vector< pair > matchFramesWithSegmentIdAndSliceNumber(FGInterface &fgInterface); + + int unpackFramesAndWriteSegmentImage(vector< pair > matchingSegmentIDsAndSliceNumbers); + int initializeMetaDataFromDICOM(DcmDataset*); + int createDICOMSegmentation(); + int createNewSegmentImage(Uint16 segmentId); + + Json::Value getSegmentAttributesMetadata(); + + IODGeneralEquipmentModule::EquipmentInfo generalEquipmentInfoModule; + + Uint16 getSegmentId(FGInterface &fgInterface, size_t frameId) const; + + template + vector > getSliceMapForSegmentation2DerivationImage(const vector dcmDatasets, + const ImageTypePointer &labelImage); + + bool hasDerivationImages(vector > &slice2derimg) const; +}; + + +#endif //DCMQI_SEGMENTATIONIMAGEOBJECT_H diff --git a/libsrc/CMakeLists.txt b/libsrc/CMakeLists.txt index 02e7caf2..29195d97 100644 --- a/libsrc/CMakeLists.txt +++ b/libsrc/CMakeLists.txt @@ -1,39 +1,16 @@ #----------------------------------------------------------------------------- -set(ADDITIONAL_SRCS) -set(INCLUDE_DIR ../include/dcmqi) - -set(HDRS - ${INCLUDE_DIR}/preproc.h - ${INCLUDE_DIR}/QIICRConstants.h - ${INCLUDE_DIR}/QIICRUIDs.h - ${INCLUDE_DIR}/ConverterBase.h - ${INCLUDE_DIR}/Exceptions.h - ${INCLUDE_DIR}/framesorter.h - ${INCLUDE_DIR}/ImageSEGConverter.h - ${INCLUDE_DIR}/ParaMapConverter - ${INCLUDE_DIR}/Helper.h - ${INCLUDE_DIR}/JSONMetaInformationHandlerBase.h - ${INCLUDE_DIR}/JSONParametricMapMetaInformationHandler.h - ${INCLUDE_DIR}/JSONSegmentationMetaInformationHandler.h - ${INCLUDE_DIR}/SegmentAttributes.h +file(GLOB_RECURSE HDRS RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} + ${PROJECT_SOURCE_DIR}/include/dcmqi/*.h ) -set(SRCS - ConverterBase.cpp - ImageSEGConverter.cpp - ParaMapConverter.cpp - Helper.cpp - JSONMetaInformationHandlerBase.cpp - JSONParametricMapMetaInformationHandler.cpp - JSONSegmentationMetaInformationHandler.cpp - SegmentAttributes.cpp +file(GLOB_RECURSE SRCS RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} + ${PROJECT_SOURCE_DIR}/libsrc/*.cpp ) - if(DCMQI_BUILTIN_JSONCPP) - list(APPEND ADDITIONAL_SRCS + list(APPEND SRCS ${DCMQI_SOURCE_DIR}/jsoncpp/jsoncpp.cpp ) set(JsonCpp_INCLUDE_DIR ${DCMQI_SOURCE_DIR}/jsoncpp) @@ -44,7 +21,6 @@ set(lib_name dcmqi) add_library(${lib_name} STATIC ${HDRS} ${SRCS} - ${ADDITIONAL_SRCS} ) if(DCMQI_LIBRARY_PROPERTIES) @@ -72,4 +48,4 @@ target_link_libraries(${lib_name} PUBLIC ${DCMTK_LIBRARIES} ${ITK_LIBRARIES} $<$>:${JsonCpp_LIBRARY}> - ) + ) \ No newline at end of file diff --git a/libsrc/DICOMFrame.cpp b/libsrc/DICOMFrame.cpp new file mode 100644 index 00000000..f652f26a --- /dev/null +++ b/libsrc/DICOMFrame.cpp @@ -0,0 +1,22 @@ +// +// Created by Andrey Fedorov on 3/9/17. +// + +#include "dcmqi/DICOMFrame.h" + +namespace dcmqi { + int DICOMFrame::initializeFrameGeometryFromLegacyInstance() { + OFString ippStr; + frameIPP.set_size(3); + for(int j=0;j<3;j++){ + CHECK_COND(frameDataset->findAndGetOFString(DCM_ImagePositionPatient, ippStr, j)); + frameIPP[j] = atof(ippStr.c_str()); + } + return EXIT_SUCCESS; + } + + int DICOMFrame::initializeFrameGeometryFromEnhancedInstance() { + OFString ippStr; + return EXIT_SUCCESS; + } +} diff --git a/libsrc/Helper.cpp b/libsrc/Helper.cpp index 92f69c3a..422db6d6 100644 --- a/libsrc/Helper.cpp +++ b/libsrc/Helper.cpp @@ -1,5 +1,6 @@ // DCMQI includes +#include #include "dcmqi/Helper.h" namespace dcmqi { @@ -368,5 +369,53 @@ namespace dcmqi { return new CodeSequenceMacro(code.c_str(), designator.c_str(), meaning.c_str()); } + OFString Helper::generateUID() { + char charUID[128]; + dcmGenerateUniqueIdentifier(charUID, QIICR_UID_ROOT); + return OFString(charUID); + } + + OFString Helper::getTagAsOFString(DcmDataset* dcm, DcmTagKey tag) { + OFString value; + CHECK_COND(dcm->findAndGetOFString(tag, value)); + return value; + } + + CodeSequenceMacro Helper::jsonToCodeSequenceMacro(Json::Value jv){ + return CodeSequenceMacro(jv["CodeValue"].asCString(), + jv["CodingSchemeDesignator"].asCString(), + jv["CodeMeaning"].asCString()); + } + + string Helper::codeSequenceMacroToString(CodeSequenceMacro c){ + OFString codeValue, codingSchemeDesignator, codeMeaning; + string s = string()+codeValue.c_str()+","+codingSchemeDesignator.c_str()+","+codeMeaning.c_str(); + return s; + } + Json::Value Helper::codeSequence2Json(CodeSequenceMacro &codeSequence) { + Json::Value value; + value["CodeValue"] = getCodeSequenceValue(codeSequence); + value["CodingSchemeDesignator"] = getCodeSequenceDesignator(codeSequence); + value["CodeMeaning"] = getCodeSequenceMeaning(codeSequence); + return value; + } + + string Helper::getCodeSequenceValue(CodeSequenceMacro &codeSequence) { + OFString value; + codeSequence.getCodeValue(value); + return value.c_str(); + } + + string Helper::getCodeSequenceDesignator(CodeSequenceMacro &codeSequence) { + OFString designator; + codeSequence.getCodingSchemeDesignator(designator); + return designator.c_str(); + } + + string Helper::getCodeSequenceMeaning(CodeSequenceMacro &codeSequence) { + OFString meaning; + codeSequence.getCodeMeaning(meaning); + return meaning.c_str(); + } } diff --git a/libsrc/ImageVolume.cpp b/libsrc/ImageVolume.cpp new file mode 100644 index 00000000..ebf64482 --- /dev/null +++ b/libsrc/ImageVolume.cpp @@ -0,0 +1,215 @@ +// +// Created by Andrey Fedorov on 3/9/17. +// + +#include +#include +#include + +#include +#include +#include +#include "dcmqi/ImageVolume.h" + +using namespace std; + +int dcmqi::ImageVolume::initializeFromDICOM(FGInterface& fgInterface) { + if(initializeDirections(fgInterface)) + return EXIT_FAILURE; + if(initializeExtent(fgInterface)) + return EXIT_FAILURE; + return EXIT_SUCCESS; +} + + +int dcmqi::ImageVolume::initializeFromITK(Float32ITKImageType::Pointer itkImage){ + return EXIT_SUCCESS; +} + +int dcmqi::ImageVolume::initializeDirections(FGInterface &fgInterface) { + // TODO: handle the situation when FoR is not initialized + OFBool isPerFrame; + + FGPlaneOrientationPatient *planorfg = OFstatic_cast(FGPlaneOrientationPatient*, + fgInterface.get(0, DcmFGTypes::EFG_PLANEORIENTPATIENT, isPerFrame)); + if(!planorfg){ + cerr << "Plane Orientation (Patient) is missing, cannot parse input " << endl; + return EXIT_FAILURE; + } + + OFString orientStr; + for(int i=0;i<3;i++){ + if(planorfg->getImageOrientationPatient(orientStr, i).good()){ + rowDirection[i] = atof(orientStr.c_str()); + } else { + cerr << "Failed to get orientation " << i << endl; + return EXIT_FAILURE; + } + } + for(int i=3;i<6;i++){ + if(planorfg->getImageOrientationPatient(orientStr, i).good()){ + columnDirection[i-3] = atof(orientStr.c_str()); + } else { + cerr << "Failed to get orientation " << i << endl; + return EXIT_FAILURE; + } + } + vnl_vector sliceDirection = vnl_cross_3d(rowDirection, columnDirection); + sliceDirection.normalize(); + + cout << "Row direction: " << rowDirection << endl; + cout << "Col direction: " << columnDirection << endl; + cout << "Z direction: " << sliceDirection << endl; + + return EXIT_SUCCESS; +} + + +int dcmqi::ImageVolume::initializeExtent(FGInterface &fgInterface) { + // Size + // Rows/Columns can be read directly from the respective attributes + // For number of slices, consider that all segments must have the same number of frames. + // If we have FoR UID initialized, this means every segment should also have Plane + // Position (Patient) initialized. So we can get the number of slices by looking + // how many per-frame functional groups a segment has. + + vector originDistances; + map originStr2distance; + map frame2overlap; + double minDistance = 0.0; + + double calculatedSliceSpacing = 0; + + unsigned numFrames = fgInterface.getNumberOfFrames(); + + /* Framesorter is to be moved to DCMTK at some point + * in the future. For now it is causing build issues on windows + + FrameSorterIPP fsIPP; + FrameSorterIPP::Results sortResults; + fsIPP.setSorterInput(&fgInterface); + fsIPP.sort(sortResults); + + */ + + // Determine ordering of the frames, keep mapping from ImagePositionPatient string + // to the distance, and keep track (just out of curiosity) how many frames overlap + vnl_vector refOrigin(3); + { + OFBool isPerFrame; + FGPlanePosPatient *planposfg = OFstatic_cast(FGPlanePosPatient*, + fgInterface.get(0, DcmFGTypes::EFG_PLANEPOSPATIENT, isPerFrame)); + if(!isPerFrame){ + cerr << "PlanePositionPatient FG is shared ... cannot handle this!" << endl; + return EXIT_FAILURE; + } + + for(int j=0;j<3;j++){ + OFString planposStr; + if(planposfg->getImagePositionPatient(planposStr, j).good()){ + refOrigin[j] = atof(planposStr.c_str()); + } else { + cerr << "Failed to read patient position" << endl; + } + } + } + + for(size_t frameId=0;frameId sOrigin; + OFString sOriginStr = ""; + sOrigin.set_size(3); + for(int j=0;j<3;j++){ + OFString planposStr; + if(planposfg->getImagePositionPatient(planposStr, j).good()){ + sOrigin[j] = atof(planposStr.c_str()); + sOriginStr += planposStr; + if(j<2) + sOriginStr+='/'; + } else { + cerr << "Failed to read patient position" << endl; + return EXIT_FAILURE; + } + } + + // check if this frame has already been encountered + if(originStr2distance.find(sOriginStr) == originStr2distance.end()){ + vnl_vector difference; + difference.set_size(3); + difference[0] = sOrigin[0]-refOrigin[0]; + difference[1] = sOrigin[1]-refOrigin[1]; + difference[2] = sOrigin[2]-refOrigin[2]; + double dist = dot_product(difference,sliceDirection); + frame2overlap[sOriginStr] = 1; + originStr2distance[sOriginStr] = dist; + assert(originStr2distance.find(sOriginStr) != originStr2distance.end()); + originDistances.push_back(dist); + + if(frameId==0){ + minDistance = dist; + origin[0] = sOrigin[0]; + origin[1] = sOrigin[1]; + origin[2] = sOrigin[2]; + } + else if(dist1){ + // WARNING: this should be improved further. Spacing should be calculated for + // consecutive frames of the individual segment. Right now, all frames are considered + // indiscriminately. Question is whether it should be computed at all, considering we do + // not have any information about whether the 2 frames are adjacent or not, so perhaps we should + // always rely on the declared spacing, and not even try to compute it? + // TODO: discuss this with the QIICR team! + + // sort all unique distances, this will be used to check consistency of + // slice spacing, and also to locate the slice position from ImagePositionPatient + // later when we read the segments + sort(originDistances.begin(), originDistances.end()); + + calculatedSliceSpacing = fabs(originDistances[0]-originDistances[1]); + + for(size_t i=1;i::const_iterator it=frame2overlap.begin(); + it!=frame2overlap.end();++it){ + if(it->second>1) + overlappingFramesCnt++; + } + + cout << "Total frames: " << numFrames << endl; + cout << "Total frames with unique IPP: " << originDistances.size() << endl; + cout << "Total overlapping frames: " << overlappingFramesCnt << endl; + cout << "Origin: " << origin << endl; + } + + return EXIT_SUCCESS; +} + diff --git a/libsrc/ImageVolumeGeometry.cpp b/libsrc/ImageVolumeGeometry.cpp new file mode 100644 index 00000000..10570b0f --- /dev/null +++ b/libsrc/ImageVolumeGeometry.cpp @@ -0,0 +1,42 @@ +// +// Created by Andrey Fedorov on 3/11/17. +// + +#include "dcmqi/ImageVolumeGeometry.h" + +ImageVolumeGeometry::ImageVolumeGeometry() { + rowDirection.set_size(3); + columnDirection.set_size(3); + sliceDirection.set_size(3); + spacing.set_size(3); + extent.set_size(3); + origin.set_size(3); +} + +int ImageVolumeGeometry::setSpacing(DoubleVectorType s) { + for(int i=0;i<3;i++) + spacing[i] = s[i]; + return EXIT_SUCCESS; +} + +int ImageVolumeGeometry::setOrigin(PointType s) { + for(int i=0;i<3;i++) + origin[i] = s[i]; + return EXIT_SUCCESS; +} + +int ImageVolumeGeometry::setExtent(SizeType s) { + for (int i = 0; i < 3; i++) + extent[i] = s[i]; + return EXIT_SUCCESS; +} + +int ImageVolumeGeometry::setDirections(DirectionType d) { + for (int i = 0; i < 3; i++) + rowDirection[i] = d[i][0]; + for (int i = 0; i < 3; i++) + columnDirection[i] = d[i][1]; + for (int i = 0; i < 3; i++) + sliceDirection[i] = d[i][2]; + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/libsrc/ConverterBase.cpp b/libsrc/MultiframeConverter.cpp similarity index 81% rename from libsrc/ConverterBase.cpp rename to libsrc/MultiframeConverter.cpp index 836209eb..1d4dce88 100644 --- a/libsrc/ConverterBase.cpp +++ b/libsrc/MultiframeConverter.cpp @@ -1,11 +1,11 @@ // DCMQI includes -#include "dcmqi/ConverterBase.h" +#include "dcmqi/MultiframeConverter.h" namespace dcmqi { - IODGeneralEquipmentModule::EquipmentInfo ConverterBase::getEquipmentInfo() { + IODGeneralEquipmentModule::EquipmentInfo MultiframeConverter::getEquipmentInfo() { // TODO: change to following for most recent dcmtk // return IODGeneralEquipmentModule::EquipmentInfo(QIICR_MANUFACTURER, QIICR_DEVICE_SERIAL_NUMBER, // QIICR_MANUFACTURER_MODEL_NAME, QIICR_SOFTWARE_VERSIONS); @@ -17,13 +17,13 @@ namespace dcmqi { return eq; } - IODEnhGeneralEquipmentModule::EquipmentInfo ConverterBase::getEnhEquipmentInfo() { + IODEnhGeneralEquipmentModule::EquipmentInfo MultiframeConverter::getEnhEquipmentInfo() { return IODEnhGeneralEquipmentModule::EquipmentInfo(QIICR_MANUFACTURER, QIICR_DEVICE_SERIAL_NUMBER, QIICR_MANUFACTURER_MODEL_NAME, QIICR_SOFTWARE_VERSIONS); } // TODO: defaults for sub classes needs to be defined - ContentIdentificationMacro ConverterBase::createContentIdentificationInformation(JSONMetaInformationHandlerBase &metaInfo) { + ContentIdentificationMacro MultiframeConverter::createContentIdentificationInformation(JSONMetaInformationHandlerBase &metaInfo) { ContentIdentificationMacro ident; CHECK_COND(ident.setContentCreatorName("dcmqi")); if(metaInfo.metaInfoRoot["seriesAttributes"].isMember("ContentDescription")){ diff --git a/libsrc/MultiframeObject.cpp b/libsrc/MultiframeObject.cpp new file mode 100644 index 00000000..20e239f3 --- /dev/null +++ b/libsrc/MultiframeObject.cpp @@ -0,0 +1,545 @@ +// +// Created by Andrey Fedorov on 3/11/17. +// + +#include +#include "dcmqi/MultiframeObject.h" + +int MultiframeObject::initializeFromDICOM(std::vector sourceDataset) { + return EXIT_SUCCESS; +} + +int MultiframeObject::initializeMetaDataFromString(const std::string &metaDataStr) { + std::istringstream metaDataStream(metaDataStr); + metaDataStream >> metaDataJson; + return EXIT_SUCCESS; +} + +int MultiframeObject::initializeContentIdentification() { + + if(sourceRepresentationType == ITK_REPR){ + CHECK_COND(contentIdentificationMacro.setContentCreatorName("dcmqi")); + if(metaDataJson.isMember("ContentDescription")){ + CHECK_COND(contentIdentificationMacro.setContentDescription(metaDataJson["ContentDescription"].asCString())); + } else { + CHECK_COND(contentIdentificationMacro.setContentDescription("DCMQI")); + } + if(metaDataJson.isMember("ContentLabel")){ + CHECK_COND(contentIdentificationMacro.setContentLabel(metaDataJson["ContentLabel"].asCString())); + } else { + CHECK_COND(contentIdentificationMacro.setContentLabel("DCMQI")); + } + if(metaDataJson.isMember("InstanceNumber")){ + CHECK_COND(contentIdentificationMacro.setInstanceNumber(metaDataJson["InstanceNumber"].asCString())); + } else { + CHECK_COND(contentIdentificationMacro.setInstanceNumber("1")); + } + CHECK_COND(contentIdentificationMacro.check()) + return EXIT_SUCCESS; + } else { // DICOM_REPR + } + return EXIT_SUCCESS; +} + +int MultiframeObject::initializeVolumeGeometryFromITK(DummyImageType::Pointer image) { + DummyImageType::SpacingType spacing; + DummyImageType::PointType origin; + DummyImageType::DirectionType directions; + DummyImageType::SizeType extent; + + spacing = image->GetSpacing(); + directions = image->GetDirection(); + extent = image->GetLargestPossibleRegion().GetSize(); + origin = image->GetOrigin(); + + volumeGeometry.setSpacing(spacing); + volumeGeometry.setOrigin(origin); + volumeGeometry.setExtent(extent); + volumeGeometry.setDirections(directions); + + return EXIT_SUCCESS; +} + +int MultiframeObject::initializeVolumeGeometryFromDICOM(FGInterface &fgInterface) { + SpacingType spacing; + PointType origin; + DirectionType directions; + SizeType extent; + + if (getImageDirections(fgInterface, directions)) { + cerr << "Failed to get image directions" << endl; + throw -1; + } + + double computedSliceSpacing, computedVolumeExtent; + vnl_vector sliceDirection(3); + sliceDirection[0] = directions[0][2]; + sliceDirection[1] = directions[1][2]; + sliceDirection[2] = directions[2][2]; + if (computeVolumeExtent(fgInterface, sliceDirection, origin, computedSliceSpacing, computedVolumeExtent)) { + cerr << "Failed to compute origin and/or slice spacing!" << endl; + throw -1; + } + + if (getDeclaredImageSpacing(fgInterface, spacing)) { + cerr << "Failed to get image spacing from DICOM!" << endl; + throw -1; + } + + const double tolerance = 1e-5; + if(!spacing[2]){ + spacing[2] = computedSliceSpacing; + } else if(fabs(spacing[2]-computedSliceSpacing)>tolerance){ + cerr << "WARNING: Declared slice spacing is significantly different from the one declared in DICOM!" << + " Declared = " << spacing[2] << " Computed = " << computedSliceSpacing << endl; + } + + OFString str; + if(dcmRepresentation->findAndGetOFString(DCM_Rows, str).good()) + extent[1] = atoi(str.c_str()); + if(dcmRepresentation->findAndGetOFString(DCM_Columns, str).good()) + extent[0] = atoi(str.c_str()); + + cout << "computed extent: " << computedVolumeExtent << "/" << spacing[2] << endl; + + // difference between the origins of the first and last frames, divide by SpacingBetweenSlices + // (which must be declared for meaningful interpretation of the SEG geometry, as amended in CP-1427), and add 1. + // see also: https://github.com/QIICR/dcmqi/issues/275 + extent[2] = round(computedVolumeExtent/spacing[2] + 1); + + cout << "Extent:" << extent << endl; + + volumeGeometry.setSpacing(spacing); + volumeGeometry.setOrigin(origin); + volumeGeometry.setExtent(extent); + volumeGeometry.setDirections(directions); + + return EXIT_SUCCESS; +} + +// for now, we do not support initialization from JSON only, +// and we don't have a way to validate metadata completeness - TODO! +bool MultiframeObject::metaDataIsComplete() { + return false; +} + +// List of tags, and FGs they belong to, for initializing dimensions module +int MultiframeObject::initializeDimensions(std::vector > dimTagList){ + OFString dimUID; + + dimensionsModule.clearData(); + + dimUID = dcmqi::Helper::generateUID(); + for(int i=0;i dimTagPair = dimTagList[i]; + CHECK_COND(dimensionsModule.addDimensionIndex(dimTagPair.first, dimUID, dimTagPair.second, + dimTagPair.first.getTagName())); + } + return EXIT_SUCCESS; +} + +int MultiframeObject::initializePixelMeasuresFG(){ + string pixelSpacingStr, sliceSpacingStr; + + pixelSpacingStr = dcmqi::Helper::floatToStrScientific(volumeGeometry.spacing[0])+ + "\\"+dcmqi::Helper::floatToStrScientific(volumeGeometry.spacing[1]); + sliceSpacingStr = dcmqi::Helper::floatToStrScientific(volumeGeometry.spacing[2]); + + CHECK_COND(pixelMeasuresFG.setPixelSpacing(pixelSpacingStr.c_str())); + CHECK_COND(pixelMeasuresFG.setSpacingBetweenSlices(sliceSpacingStr.c_str())); + CHECK_COND(pixelMeasuresFG.setSliceThickness(sliceSpacingStr.c_str())); + + return EXIT_SUCCESS; +} + +int MultiframeObject::initializePlaneOrientationFG() { + planeOrientationPatientFG.setImageOrientationPatient( + dcmqi::Helper::floatToStrScientific(volumeGeometry.rowDirection[0]).c_str(), + dcmqi::Helper::floatToStrScientific(volumeGeometry.rowDirection[1]).c_str(), + dcmqi::Helper::floatToStrScientific(volumeGeometry.rowDirection[2]).c_str(), + dcmqi::Helper::floatToStrScientific(volumeGeometry.columnDirection[0]).c_str(), + dcmqi::Helper::floatToStrScientific(volumeGeometry.columnDirection[1]).c_str(), + dcmqi::Helper::floatToStrScientific(volumeGeometry.columnDirection[2]).c_str() + ); + return EXIT_SUCCESS; +} + +ContentItemMacro* MultiframeObject::initializeContentItemMacro(CodeSequenceMacro conceptName, + CodeSequenceMacro conceptCode){ + ContentItemMacro* item = new ContentItemMacro(); + CodeSequenceMacro* concept = new CodeSequenceMacro(conceptName); + CodeSequenceMacro* value = new CodeSequenceMacro(conceptCode); + + if (!item || !concept || !value) + { + return NULL; + } + + item->getEntireConceptNameCodeSequence().push_back(concept); + item->getEntireConceptCodeSequence().push_back(value); + item->setValueType(ContentItemMacro::VT_CODE); + + return EXIT_SUCCESS; +} + +// populates slice2frame vector that maps each of the volume slices to the set of frames that +// are considered as derivation dataset +// TODO: this function assumes that all of the derivation datasets are images, which is probably ok +int MultiframeObject::mapVolumeSlicesToDICOMFrames(const vector dcmDatasets, + vector > &slice2frame){ + slice2frame.resize(volumeGeometry.extent[2]); + + for(vector::const_iterator it=dcmDatasets.begin();it!=dcmDatasets.end();++it) { + Uint32 numFrames; + DcmDataset* dcm = *it; + if(dcm->findAndGetUint32(DCM_NumberOfFrames, numFrames).good()){ + // this is a multiframe object + for(int f=0;f > &slice2frame, DcmDataset *dcm, int frameNo) { + dcmqi::DICOMFrame frame(dcm, frameNo); + vector intersectingSlices = findIntersectingSlices(volumeGeometry, frame); + + if(intersectingSlices.size()) + insertDerivationSeriesInstance(frame.getSeriesUID(), frame.getInstanceUID()); + + for(int s=0;s derivationFrames, + CodeSequenceMacro purposeOfReferenceCode, + CodeSequenceMacro derivationCode){ + DerivationImageItem *derimgItem; + // TODO: check what we should do with DerivationDescription + CHECK_COND(fgder->addDerivationImageItem(derivationCode,"",derimgItem)); + + OFVector siVector; + std::vector frameVector; + + for(set::const_iterator sIt=derivationFrames.begin(); + sIt!=derivationFrames.end();++sIt) { + siVector.push_back((*sIt).getDataset()); + frameVector.push_back(*sIt); + } + + OFVector srcimgItems; + CHECK_COND(derimgItem->addSourceImageItems(siVector, purposeOfReferenceCode, srcimgItems)); + + // iterate over source image items (assuming they are in the same order as in siVector!), and initialize + // frame number, if applicable + unsigned siItemCnt=0; + for(OFVector::iterator vIt=srcimgItems.begin(); + vIt!=srcimgItems.end();++vIt,++siItemCnt) { + // TODO: when multuple frames from the same instance are used, they should be referenced within a single + // ImageSOPInstanceReferenceMacro. There would need to be another level of checks over all of the frames + // that are mapped to the given slice to identify those that are from the same instance, and populate the + // list of frames. I can't think of any use case where this would be immediately important, but if we ever + // use multiframe for DCE/DWI, with all of the temporal/b-value frames stored in a single object, there will + // be multiple frames used to derive a single frame in a parametric map, for example. + if(frameVector[siItemCnt].getFrameNumber()) { + OFVector frameNumbersVector; + ImageSOPInstanceReferenceMacro &instRef = (*vIt)->getImageSOPInstanceReference(); + frameNumbersVector.push_back(frameVector[siItemCnt].getFrameNumber()); + instRef.setReferencedFrameNumber(frameNumbersVector); + } + } + + return EXIT_SUCCESS; +} + +std::vector MultiframeObject::findIntersectingSlices(ImageVolumeGeometry &volume, dcmqi::DICOMFrame &frame) { + std::vector intersectingSlices; + // for now, adopt a simple strategy that maps origin of the frame to index, and selects the slice corresponding + // to this index as the intersecting one + ImageVolumeGeometry::DummyImageType::Pointer itkvolume = volume.getITKRepresentation(); + ImageVolumeGeometry::DummyImageType::PointType point; + ImageVolumeGeometry::DummyImageType::IndexType index; + vnl_vector frameIPP = frame.getFrameIPP(); + point[0] = frameIPP[0]; + point[1] = frameIPP[1]; + point[2] = frameIPP[2]; + + if(itkvolume->TransformPhysicalPointToIndex(point, index)) { + intersectingSlices.push_back(index[2]); + } + + return intersectingSlices; +} + +void MultiframeObject::insertDerivationSeriesInstance(string seriesUID, string instanceUID){ + if(derivationSeriesToInstanceUIDs.find(seriesUID) == derivationSeriesToInstanceUIDs.end()) + derivationSeriesToInstanceUIDs[seriesUID] = std::set(); + derivationSeriesToInstanceUIDs[seriesUID].insert(instanceUID); +} + +int MultiframeObject::initializeCommonInstanceReferenceModule(IODCommonInstanceReferenceModule &commref, vector > &slice2frame){ + + // map individual series UIDs to the list of instance UIDs - we need to have this organization + // to populate Common instance reference module + std::map > series2frame; + for(int slice=0;slice!=slice2frame.size();slice++){ + for(set::const_iterator frameI=slice2frame[slice].begin(); + frameI!=slice2frame[slice].end();++frameI){ + dcmqi::DICOMFrame frame = *frameI; + if(series2frame.find(frame.getSeriesUID()) == series2frame.end()){ + std::set setOfInstances; + setOfInstances.insert(frame); + series2frame[frame.getSeriesUID()] = setOfInstances; + } else { + series2frame[frame.getSeriesUID()].insert(frame); + } + } + } + + // create a new ReferencedSeriesItem for each series, and populate with instances + OFVector &refseries = commref.getReferencedSeriesItems(); + for(std::map >::const_iterator mIt=series2frame.begin(); mIt!=series2frame.end();++mIt){ + // TODO: who is supposed to de-allocate this? + IODSeriesAndInstanceReferenceMacro::ReferencedSeriesItem* refseriesItem = new IODSeriesAndInstanceReferenceMacro::ReferencedSeriesItem; + refseriesItem->setSeriesInstanceUID(mIt->first.c_str()); + OFVector &refinstances = refseriesItem->getReferencedInstanceItems(); + + for(std::set::const_iterator sIt=mIt->second.begin();sIt!=mIt->second.end();++sIt){ + dcmqi::DICOMFrame frame = *sIt; + SOPInstanceReferenceMacro *refinstancesItem = new SOPInstanceReferenceMacro(); + CHECK_COND(refinstancesItem->setReferencedSOPClassUID(frame.getClassUID().c_str())); + CHECK_COND(refinstancesItem->setReferencedSOPInstanceUID(frame.getInstanceUID().c_str())); + refinstances.push_back(refinstancesItem); + } + + refseries.push_back(refseriesItem); + } + + return EXIT_SUCCESS; +} + +int MultiframeObject::getImageDirections(FGInterface& fgInterface, DirectionType &dir){ + // TODO: handle the situation when FoR is not initialized + OFBool isPerFrame; + vnl_vector rowDirection(3), colDirection(3); + + FGPlaneOrientationPatient *planorfg = OFstatic_cast(FGPlaneOrientationPatient*, + fgInterface.get(0, DcmFGTypes::EFG_PLANEORIENTPATIENT, isPerFrame)); + if(!planorfg){ + cerr << "Plane Orientation (Patient) is missing, cannot parse input " << endl; + return EXIT_FAILURE; + } + + if(planorfg->getImageOrientationPatient(rowDirection[0], rowDirection[1], rowDirection[2], + colDirection[0], colDirection[1], colDirection[2]).bad()){ + cerr << "Failed to get orientation " << endl; + return EXIT_FAILURE; + } + + vnl_vector sliceDirection = vnl_cross_3d(rowDirection, colDirection); + sliceDirection.normalize(); + + for(int i=0;i<3;i++){ + dir[i][0] = rowDirection[i]; + dir[i][1] = colDirection[i]; + dir[i][2] = sliceDirection[i]; + } + + cout << "Row direction: " << rowDirection << endl; + cout << "Col direction: " << colDirection << endl; + cout << "Z direction : " << sliceDirection << endl; + + return EXIT_SUCCESS; +} + +int MultiframeObject::computeVolumeExtent(FGInterface& fgInterface, vnl_vector &sliceDirection, + PointType &imageOrigin, double &sliceSpacing, double &volumeExtent) { + // Size + // Rows/Columns can be read directly from the respective attributes + // For number of slices, consider that all segments must have the same number of frames. + // If we have FoR UID initialized, this means every segment should also have Plane + // Position (Patient) initialized. So we can get the number of slices by looking + // how many per-frame functional groups a segment has. + + vector originDistances; + map originStr2distance; + map frame2overlap; + double minDistance = 0.0; + + sliceSpacing = 0; + + // Determine ordering of the frames, keep mapping from ImagePositionPatient string + // to the distance, and keep track (just out of curiosity) how many frames overlap + vnl_vector refOrigin = getFrameOrigin(fgInterface, 0); + + unsigned numFrames = fgInterface.getNumberOfFrames(); + + for(size_t frameId=0;frameId sOrigin(3); + sOrigin = getFrameOrigin(planposfg); + std::ostringstream sstream; + sstream << sOrigin[0] << "/" << sOrigin[1] << "/" << sOrigin[2]; + OFString sOriginStr = sstream.str().c_str(); + + // check if this frame has already been encountered + if(originStr2distance.find(sOriginStr) == originStr2distance.end()){ + vnl_vector difference(3); + difference[0] = sOrigin[0]-refOrigin[0]; + difference[1] = sOrigin[1]-refOrigin[1]; + difference[2] = sOrigin[2]-refOrigin[2]; + + double dist = dot_product(difference, sliceDirection); + frame2overlap[sOriginStr] = 1; + originStr2distance[sOriginStr] = dist; + assert(originStr2distance.find(sOriginStr) != originStr2distance.end()); + originDistances.push_back(dist); + + if(frameId == 0 || dist < minDistance){ + minDistance = dist; + imageOrigin[0] = sOrigin[0]; + imageOrigin[1] = sOrigin[1]; + imageOrigin[2] = sOrigin[2]; + } + } else { + frame2overlap[sOriginStr]++; + } + } + + // it IS possible to have a segmentation object containing just one frame! + if(numFrames>1){ + // WARNING: this should be improved further. Spacing should be calculated for + // consecutive frames of the individual segment. Right now, all frames are considered + // indiscriminately. Question is whether it should be computed at all, considering we do + // not have any information about whether the 2 frames are adjacent or not, so perhaps we should + // always rely on the declared spacing, and not even try to compute it? + // TODO: discuss this with the QIICR team! + + // sort all unique distances, this will be used to check consistency of + // slice spacing, and also to locate the slice position from ImagePositionPatient + // later when we read the segments + sort(originDistances.begin(), originDistances.end()); + + sliceSpacing = fabs(originDistances[0]-originDistances[1]); + + for(size_t i=1;i::const_iterator it=frame2overlap.begin(); + it!=frame2overlap.end();++it){ + if(it->second>1) + overlappingFramesCnt++; + } + + cout << "Total frames: " << numFrames << endl; + cout << "Total frames with unique IPP: " << originDistances.size() << endl; + cout << "Total overlapping frames: " << overlappingFramesCnt << endl; + cout << "Origin: " << imageOrigin << endl; + } + + return EXIT_SUCCESS; +} + +vnl_vector MultiframeObject::getFrameOrigin(FGInterface &fgInterface, int frameId) const { + OFBool isPerFrame; + FGPlanePosPatient *planposfg = OFstatic_cast(FGPlanePosPatient*, + fgInterface.get(frameId, DcmFGTypes::EFG_PLANEPOSPATIENT, isPerFrame)); + return getFrameOrigin(planposfg); +} + +vnl_vector MultiframeObject::getFrameOrigin(FGPlanePosPatient *planposfg) const { + vnl_vector origin(3); + if(!planposfg->getImagePositionPatient(origin[0], origin[1], origin[2]).good()){ + cerr << "Failed to read patient position" << endl; + throw -1; + } + return origin; +} + +int MultiframeObject::getDeclaredImageSpacing(FGInterface &fgInterface, SpacingType &spacing){ + OFBool isPerFrame; + FGPixelMeasures *pixelMeasures = OFstatic_cast(FGPixelMeasures*, + fgInterface.get(0, DcmFGTypes::EFG_PIXELMEASURES, isPerFrame)); + if(!pixelMeasures){ + cerr << "Pixel measures FG is missing!" << endl; + return EXIT_FAILURE; + } + + pixelMeasures->getPixelSpacing(spacing[0], 0); + pixelMeasures->getPixelSpacing(spacing[1], 1); + + Float64 spacingFloat; + if(pixelMeasures->getSpacingBetweenSlices(spacingFloat,0).good() && spacingFloat != 0){ + spacing[2] = spacingFloat; + } else if(pixelMeasures->getSliceThickness(spacingFloat,0).good() && spacingFloat != 0){ + // SliceThickness can be carried forward from the source images, and may not be what we need + // As an example, this ePAD example has 1.25 carried from CT, but true computed thickness is 1! + cerr << "WARNING: SliceThickness is present and is " << spacingFloat << ". using it!" << endl; + spacing[2] = spacingFloat; + } + return EXIT_SUCCESS; +} + +int MultiframeObject::initializeSeriesSpecificAttributes(IODGeneralSeriesModule& generalSeriesModule, + IODGeneralImageModule& generalImageModule){ + + // TODO: error checks + string bodyPartExamined; + if(metaDataJson.isMember("BodyPartExamined")){ + bodyPartExamined = metaDataJson["BodyPartExamined"].asCString(); + } + if(derivationDcmDatasets.size() && bodyPartExamined.empty()) { + OFString bodyPartStr; + DcmDataset *srcDataset = derivationDcmDatasets[0]; + if(srcDataset->findAndGetOFString(DCM_BodyPartExamined, bodyPartStr).good()) { + if (!bodyPartStr.empty()) + bodyPartExamined = bodyPartStr.c_str(); + } + } + + if(!bodyPartExamined.empty()) + generalSeriesModule.setBodyPartExamined(bodyPartExamined.c_str()); + + // SeriesDate/Time should be of when parametric map was taken; initialize to when it was saved + { + OFString contentDate, contentTime; + DcmDate::getCurrentDate(contentDate); + DcmTime::getCurrentTime(contentTime); + + // TODO: AcquisitionTime + generalSeriesModule.setSeriesDate(contentDate.c_str()); + generalSeriesModule.setSeriesTime(contentTime.c_str()); + generalImageModule.setContentDate(contentDate.c_str()); + generalImageModule.setContentTime(contentTime.c_str()); + } + + generalSeriesModule.setSeriesDescription(metaDataJson["SeriesDescription"].asCString()); + generalSeriesModule.setSeriesNumber(metaDataJson["SeriesNumber"].asCString()); + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/libsrc/ParaMapConverter.cpp b/libsrc/ParametricMapConverter.cpp similarity index 90% rename from libsrc/ParaMapConverter.cpp rename to libsrc/ParametricMapConverter.cpp index 4f674d00..88379262 100644 --- a/libsrc/ParaMapConverter.cpp +++ b/libsrc/ParametricMapConverter.cpp @@ -4,54 +4,75 @@ #include // DCMQI includes -#include "dcmqi/ParaMapConverter.h" -#include "dcmqi/ImageSEGConverter.h" +#include "dcmqi/ParametricMapConverter.h" +#include "dcmqi/SegmentationImageConverter.h" +#include "dcmqi/ParametricMapObject.h" using namespace std; namespace dcmqi { - DcmDataset* ParaMapConverter::itkimage2paramap(const FloatImageType::Pointer ¶metricMapImage, vector dcmDatasets, + DcmDataset* itkimage2paramapReplacement(const FloatImageType::Pointer ¶metricMapImage, vector dcmDatasets, const string &metaData) { + ParametricMapObject pm; + pm.initializeFromITK(parametricMapImage, metaData, dcmDatasets); - MinMaxCalculatorType::Pointer calculator = MinMaxCalculatorType::New(); - calculator->SetImage(parametricMapImage); - calculator->Compute(); + DcmDataset *output = new DcmDataset(); + pm.getDICOMRepresentation(*output); + + return output; + } + + pair paramap2itkimageReplacement(DcmDataset *pmapDataset){ + + ParametricMapObject pm; + pm.initializeFromDICOM(pmapDataset); + + Json::StyledWriter styledWriter; + std::stringstream ss; + + ss << styledWriter.write(pm.getMetaDataJson()); + + return pair (pm.getITKRepresentation(), + ss.str()); + }; + + DcmDataset* ParametricMapConverter::itkimage2paramap(const FloatImageType::Pointer ¶metricMapImage, vector dcmDatasets, + const string &metaData) { JSONParametricMapMetaInformationHandler metaInfo(metaData); metaInfo.read(); - metaInfo.setFirstValueMapped(calculator->GetMinimum()); - metaInfo.setLastValueMapped(calculator->GetMaximum()); - - IODEnhGeneralEquipmentModule::EquipmentInfo eq = getEnhEquipmentInfo(); - ContentIdentificationMacro contentID = createContentIdentificationInformation(metaInfo); + // Prepare ContentIdentification + IODEnhGeneralEquipmentModule::EquipmentInfo eq = ParametricMapConverter::getEnhEquipmentInfo(); + ContentIdentificationMacro contentID = ParametricMapConverter::createContentIdentificationInformation(metaInfo); CHECK_COND(contentID.setInstanceNumber(metaInfo.getInstanceNumber().c_str())); - // TODO: following should maybe be moved to meta info - OFString imageFlavor = "VOLUME"; - OFString pixContrast = "NONE"; - if(metaInfo.metaInfoRoot.isMember("DerivedPixelContrast")){ - pixContrast = metaInfo.metaInfoRoot["DerivedPixelContrast"].asCString(); + // Get image modality from the source dataset + DcmDataset* srcDataset = NULL; + if(dcmDatasets.size()) { + srcDataset = dcmDatasets[0]; + } else { + return NULL; } + OFString modality; + CHECK_COND(srcDataset->findAndGetOFString(DCM_Modality, modality)); - // TODO: initialize modality from the source / add to schema? - OFString modality = "MR"; - + // Get size of the image to be converted FloatImageType::SizeType inputSize = parametricMapImage->GetBufferedRegion().GetSize(); cout << "Input image size: " << inputSize << endl; + // create Parametric map object OFvariant obj = DPMParametricMapIOD::create(modality, metaInfo.getSeriesNumber().c_str(), metaInfo.getInstanceNumber().c_str(), inputSize[1], inputSize[0], eq, contentID, - imageFlavor, pixContrast, DPMTypes::CQ_RESEARCH); + "VOLUME", "QUANTITY", DPMTypes::CQ_RESEARCH); if (OFCondition* pCondition = OFget(&obj)) return NULL; DPMParametricMapIOD* pMapDoc = OFget(&obj); - DcmDataset* srcDataset = NULL; if(dcmDatasets.size()){ srcDataset = dcmDatasets[0]; } @@ -61,12 +82,13 @@ namespace dcmqi { CHECK_COND(pMapDoc->import(*srcDataset, OFTrue, OFTrue, OFTrue, OFFalse)); /* Initialize dimension module */ - char dimUID[128]; - dcmGenerateUniqueIdentifier(dimUID, QIICR_UID_ROOT); IODMultiframeDimensionModule &mfdim = pMapDoc->getIODMultiframeDimensionModule(); - OFCondition result = mfdim.addDimensionIndex(DCM_ImagePositionPatient, dimUID, + // TODO: this is probably incorrect + OFCondition result = mfdim.addDimensionIndex(DCM_ImagePositionPatient, Helper::generateUID(), DCM_RealWorldValueMappingSequence, "Frame position"); + // Initialize PM geometry + // Shared FGs: PixelMeasuresSequence { FGPixelMeasures *pixmsr = new FGPixelMeasures(); @@ -80,6 +102,8 @@ namespace dcmqi { spacingSStream << scientific << labelSpacing[2]; CHECK_COND(pixmsr->setSpacingBetweenSlices(spacingSStream.str().c_str())); CHECK_COND(pixmsr->setSliceThickness(spacingSStream.str().c_str())); + + // SHARED CHECK_COND(pMapDoc->addForAllFrames(*pixmsr)); } @@ -101,9 +125,12 @@ namespace dcmqi { Helper::floatToStrScientific(labelDirMatrix[2][1]).c_str()); //CHECK_COND(planor->setImageOrientationPatient(imageOrientationPatientStr)); + // SHARED CHECK_COND(pMapDoc->addForAllFrames(*planor)); } + + // Initialize metadata - anatomy FGFrameAnatomy frameAnaFG; frameAnaFG.setLaterality(FGFrameAnatomy::str2Laterality(metaInfo.getFrameLaterality().c_str())); if(metaInfo.metaInfoRoot.isMember("AnatomicRegionSequence")){ @@ -114,6 +141,8 @@ namespace dcmqi { } else { frameAnaFG.getAnatomy().getAnatomicRegion().set("T-D0050", "SRT", "Tissue"); } + + // SHARED CHECK_COND(pMapDoc->addForAllFrames(frameAnaFG)); FGIdentityPixelValueTransformation idTransFG; @@ -123,8 +152,11 @@ namespace dcmqi { FGParametricMapFrameType frameTypeFG; std::string frameTypeStr = "DERIVED\\PRIMARY\\VOLUME\\QUANTITY"; frameTypeFG.setFrameType(frameTypeStr.c_str()); + + // SHARED CHECK_COND(pMapDoc->addForAllFrames(frameTypeFG)); + // Initialize RWVM FGRealWorldValueMapping rwvmFG; FGRealWorldValueMapping::RWVMItem* realWorldValueMappingItem = new FGRealWorldValueMapping::RWVMItem(); if (!realWorldValueMappingItem ) @@ -135,8 +167,13 @@ namespace dcmqi { realWorldValueMappingItem->setRealWorldValueSlope(metaInfo.getRealWorldValueSlope()); realWorldValueMappingItem->setRealWorldValueIntercept(atof(metaInfo.getRealWorldValueIntercept().c_str())); - realWorldValueMappingItem->setRealWorldValueFirstValueMappedSigned(metaInfo.getFirstValueMapped()); - realWorldValueMappingItem->setRealWorldValueLastValueMappedSigned(metaInfo.getLastValueMapped()); + // Calculate intensity range - required + MinMaxCalculatorType::Pointer calculator = MinMaxCalculatorType::New(); + calculator->SetImage(parametricMapImage); + calculator->Compute(); + + realWorldValueMappingItem->setRealWorldValueFirstValueMappedSigned(calculator->GetMinimum()); + realWorldValueMappingItem->setRealWorldValueLastValueMappedSigned(calculator->GetMaximum()); CodeSequenceMacro* measurementUnitCode = metaInfo.getMeasurementUnitsCode(); if (measurementUnitCode != NULL) { @@ -148,6 +185,7 @@ namespace dcmqi { // TODO: LutExplanation and LUTLabel should be added as Metainformation realWorldValueMappingItem->setLUTExplanation(metaInfo.metaInfoRoot["MeasurementUnitsCode"]["CodeMeaning"].asCString()); realWorldValueMappingItem->setLUTLabel(metaInfo.metaInfoRoot["MeasurementUnitsCode"]["CodeValue"].asCString()); + ContentItemMacro* quantity = new ContentItemMacro; CodeSequenceMacro* qCodeName = new CodeSequenceMacro("G-C1C6", "SRT", "Quantity"); CodeSequenceMacro* qSpec = new CodeSequenceMacro( @@ -165,7 +203,7 @@ namespace dcmqi { realWorldValueMappingItem->getEntireQuantityDefinitionSequence().push_back(quantity); quantity->setValueType(ContentItemMacro::VT_CODE); - // initialize optional items, if available + // initialize optional RWVM items, if available if(metaInfo.metaInfoRoot.isMember("MeasurementMethodCode")){ ContentItemMacro* measureMethod = new ContentItemMacro; CodeSequenceMacro* qCodeName = new CodeSequenceMacro("G-C306", "SRT", "Measurement Method"); @@ -187,7 +225,7 @@ namespace dcmqi { if(metaInfo.metaInfoRoot.isMember("ModelFittingMethodCode")){ ContentItemMacro* fittingMethod = new ContentItemMacro; - CodeSequenceMacro* qCodeName = new CodeSequenceMacro("DWMPxxxxx2", "99QIICR", "Model fitting method"); + CodeSequenceMacro* qCodeName = new CodeSequenceMacro("113241", "DCM", "Model fitting method"); CodeSequenceMacro* qSpec = new CodeSequenceMacro( metaInfo.metaInfoRoot["ModelFittingMethodCode"]["CodeValue"].asCString(), metaInfo.metaInfoRoot["ModelFittingMethodCode"]["CodingSchemeDesignator"].asCString(), @@ -208,7 +246,7 @@ namespace dcmqi { for(int bvalId=0;bvalIdGetOutput()); cout << "Mapping from the ITK image slices to the DICOM instances in the input list" << endl; for(int i=0;isetReferencedSOPClassUID(classUID)); @@ -375,16 +416,10 @@ namespace dcmqi { // Frame Content OFCondition result = fgfc->setDimensionIndexValues(sliceNumber+1 /* value within dimension */, 0 /* first dimension */); -#if ADD_DERIMG - // Already pushed above if siVector.size > 0 - // if(fgder) - // perFrameFGs.push_back(fgder); -#endif - DPMParametricMapIOD::FramesType frames = pMapDoc->getFrames(); result = OFget >(&frames)->addFrame(&*data.begin(), frameSize, perFrameFGs); - cout << "Frame " << sliceNumber << " added" << endl; +// cout << "Frame " << sliceNumber << " added" << endl; } // remove derivation image FG from the per-frame FGs, only if applicable! @@ -434,7 +469,7 @@ namespace dcmqi { return output; } - pair ParaMapConverter::paramap2itkimage(DcmDataset *pmapDataset) { + pair ParametricMapConverter::paramap2itkimage(DcmDataset *pmapDataset) { DcmRLEDecoderRegistration::registerCodecs(); @@ -550,7 +585,9 @@ namespace dcmqi { return pair (pmImage, metaInfo.getJSONOutputAsString()); } - OFCondition ParaMapConverter::addFrame(DPMParametricMapIOD &map, const FloatImageType::Pointer ¶metricMapImage, + // appears to be not used +#if 0 + OFCondition ParametricMapConverter::addFrame(DPMParametricMapIOD &map, const FloatImageType::Pointer ¶metricMapImage, const JSONParametricMapMetaInformationHandler &metaInfo, const unsigned long frameNo, OFVector groups) { @@ -609,8 +646,9 @@ namespace dcmqi { } return result; } +#endif - void ParaMapConverter::populateMetaInformationFromDICOM(DcmDataset *pmapDataset, + void ParametricMapConverter::populateMetaInformationFromDICOM(DcmDataset *pmapDataset, JSONParametricMapMetaInformationHandler &metaInfo) { OFvariant result = DPMParametricMapIOD::loadDataset(*pmapDataset); @@ -647,7 +685,7 @@ namespace dcmqi { Float64 slope; // TODO: replace the following call by following getter once it is available -// item->getRealWorldValueSlope(slope); + //item->getRealWorldValueSlope(slope); item->getData().findAndGetFloat64(DCM_RealWorldValueSlope, slope); metaInfo.setRealWorldValueSlope(slope); @@ -684,7 +722,7 @@ namespace dcmqi { } FGDerivationImage* derivationImage = OFstatic_cast(FGDerivationImage*, fg.get(0, DcmFGTypes::EFG_DERIVATIONIMAGE)); - + if(derivationImage){ OFVector& derivationImageItems = derivationImage->getDerivationImageItems(); if(derivationImageItems.size()>0){ diff --git a/libsrc/ParametricMapObject.cpp b/libsrc/ParametricMapObject.cpp new file mode 100644 index 00000000..cee593ef --- /dev/null +++ b/libsrc/ParametricMapObject.cpp @@ -0,0 +1,584 @@ +// +// Created by Andrey Fedorov on 3/11/17. +// + +#include +#include "dcmqi/QIICRConstants.h" +#include "dcmqi/ParametricMapObject.h" + +int ParametricMapObject::initializeFromITK(Float32ITKImageType::Pointer inputImage, + const string &metaDataStr, + std::vector derivationDatasets) { + + setDerivationDatasets(derivationDatasets); + + sourceRepresentationType = ITK_REPR; + + itkImage = inputImage; + + initializeMetaDataFromString(metaDataStr); + + if(!metaDataIsComplete()){ + updateMetaDataFromDICOM(derivationDatasets); + } + + initializeVolumeGeometry(); + + // NB: the sequence of steps initializing different components of the object parallels that + // in the original converter function. It probably makes sense to revisit the sequence + // of these steps. It does not necessarily need to happen in this order. + initializeEquipmentInfo(); + initializeContentIdentification(); + + // TODO: consider creating parametric map object after all FGs are initialized instead + createDICOMParametricMap(); + + // populate metadata about patient/study, from derivation + // datasets or from metadata + initializeCompositeContext(); + initializeSeriesSpecificAttributes(parametricMap->getIODGeneralSeriesModule(), + parametricMap->getIODGeneralImageModule()); + + // populate functional groups + std::vector > dimensionTags; + dimensionTags.push_back(std::pair(DCM_ImagePositionPatient, DCM_PlanePositionSequence)); + initializeDimensions(dimensionTags); + + initializePixelMeasuresFG(); + CHECK_COND(parametricMap->addForAllFrames(pixelMeasuresFG)); + + initializePlaneOrientationFG(); + CHECK_COND(parametricMap->addForAllFrames(planeOrientationPatientFG)); + + // PM-specific FGs + initializeFrameAnatomyFG(); + initializeRWVMFG(); + + // Mapping from parametric map volume slices to the DICOM frames + vector > slice2frame; + + // initialize referenced instances + // this is done using this utility function from the parent class, since this functionality will + // be needed both in the PM and SEG objects + mapVolumeSlicesToDICOMFrames(derivationDatasets, slice2frame); + + initializeCommonInstanceReferenceModule(this->parametricMap->getCommonInstanceReference(), slice2frame); + + initializeFrames(slice2frame); + + return EXIT_SUCCESS; +} + +int ParametricMapObject::initializeVolumeGeometry() { + if(sourceRepresentationType == ITK_REPR){ + + Float32ToDummyCasterType::Pointer caster = + Float32ToDummyCasterType::New(); + caster->SetInput(itkImage); + caster->Update(); + + MultiframeObject::initializeVolumeGeometryFromITK(caster->GetOutput()); + } else { + MultiframeObject::initializeVolumeGeometryFromDICOM(parametricMap->getFunctionalGroups()); + } + return EXIT_SUCCESS; +} + +int ParametricMapObject::updateMetaDataFromDICOM(std::vector dcmList) { + if(!dcmList.size()) + return EXIT_FAILURE; + + DcmDataset* dcm = dcmList[0]; + metaDataJson["Modality"] = + std::string(dcmqi::Helper::getTagAsOFString(dcm, DCM_Modality).c_str()); + + return EXIT_SUCCESS; +} + +int ParametricMapObject::initializeEquipmentInfo() { + if(sourceRepresentationType == ITK_REPR){ + enhancedEquipmentInfoModule = IODEnhGeneralEquipmentModule::EquipmentInfo(QIICR_MANUFACTURER, + QIICR_DEVICE_SERIAL_NUMBER, + QIICR_MANUFACTURER_MODEL_NAME, + QIICR_SOFTWARE_VERSIONS); + /* + enhancedEquipmentInfoModule.m_Manufacturer = QIICR_MANUFACTURER; + enhancedEquipmentInfoModule.m_DeviceSerialNumber = QIICR_DEVICE_SERIAL_NUMBER; + enhancedEquipmentInfoModule.m_ManufacturerModelName = QIICR_MANUFACTURER_MODEL_NAME; + enhancedEquipmentInfoModule.m_SoftwareVersions = QIICR_SOFTWARE_VERSIONS; + */ + + } else { // DICOM_REPR + } + return EXIT_SUCCESS; +} + +int ParametricMapObject::createDICOMParametricMap() { + + // create Parametric map object + + // TODO: revisit intialization of the modality - if source images are available, modality should match + // that in the source images! + OFvariant obj = + DPMParametricMapIOD::create(metaDataJson["Modality"].asCString(), + metaDataJson["SeriesNumber"].asCString(), + metaDataJson["InstanceNumber"].asCString(), + volumeGeometry.extent[1], + volumeGeometry.extent[0], + enhancedEquipmentInfoModule, + contentIdentificationMacro, + "VOLUME", "QUANTITY", + DPMTypes::CQ_RESEARCH); + // TODO: look into the following, check with @che85 on the purpose of this line! + if (OFCondition* pCondition = OFget(&obj)) { + std::cerr << "Failed to create parametric map object!" << std::endl; + return EXIT_FAILURE; + } + + parametricMap = new DPMParametricMapIOD( *OFget(&obj) ); + + // These FG are constant + FGIdentityPixelValueTransformation idTransFG; + CHECK_COND(parametricMap->addForAllFrames(idTransFG)); + + FGParametricMapFrameType frameTypeFG; + std::string frameTypeStr = "DERIVED\\PRIMARY\\VOLUME\\QUANTITY"; + frameTypeFG.setFrameType(frameTypeStr.c_str()); + CHECK_COND(parametricMap->addForAllFrames(frameTypeFG)); + + /* Initialize dimension module */ + IODMultiframeDimensionModule &mfdim = parametricMap->getIODMultiframeDimensionModule(); + OFCondition result = mfdim.addDimensionIndex(DCM_ImagePositionPatient, dcmqi::Helper::generateUID(), + DCM_PlanePositionSequence, "ImagePositionPatient"); + + return EXIT_SUCCESS; +} + +int ParametricMapObject::createITKParametricMap() { + + initializeVolumeGeometry(); + + // Initialize the image + itkImage = volumeGeometry.getITKRepresentation(); + + DPMParametricMapIOD::FramesType obj = parametricMap->getFrames(); + if (OFCondition* pCondition = OFget(&obj)) { + throw -1; + } + + FGInterface &fgInterface = parametricMap->getFunctionalGroups(); + DPMParametricMapIOD::Frames frames = *OFget >(&obj); + + createITKImageFromFrames(fgInterface, frames); + return EXIT_SUCCESS; +} + +int ParametricMapObject::createITKImageFromFrames(FGInterface &fgInterface, + DPMParametricMapIOD::Frames frames) { + for(int frameId=0;frameIdSetPixel(index, frame[pixelPosition]); + } + } + } + return EXIT_SUCCESS; +} + +int ParametricMapObject::initializeCompositeContext() { + // TODO: should this be done in the parent? + if(derivationDcmDatasets.size()){ + /* Import GeneralSeriesModule - content for the reference + from include/dcmtk/dcmiod/modgeneralseries.h: + * Modality: (CS, 1, 1) + * Series Instance Number: (UI, 1, 1) + * Series Number: (IS, 1, 2) + * Laterality: (CS, 1, 2C) + * Series Date: (DA, 1, 3) + * Series Time: (TM, 1, 3) + * Performing Physician's Name: (PN, 1, 3) + * Protocol Name: (LO, 1, 3) + * Series Description: (LO, 1, 3) + * Operators' Name: (PN, 1-n, 3) + * Body Part Examined: (CS, 1, 3) + * Patient Position: (CS, 1, 2C) + */ + CHECK_COND(parametricMap->import(*derivationDcmDatasets[0], + OFTrue, // Patient + OFTrue, // Study + OFTrue, // Frame of reference + OFFalse)); // Series + + } else { + // TODO: once we support passing of composite context in metadata, propagate it + // into parametricMap here + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} + +int ParametricMapObject::initializeFrameAnatomyFG() { + if(metaDataJson.isMember("FrameLaterality")) + frameAnatomyFG.setLaterality(FGFrameAnatomy::str2Laterality(metaDataJson["FrameLaterality"].asCString())); + else + frameAnatomyFG.setLaterality(FGFrameAnatomy::str2Laterality("U")); + + // TODO: simplify code initialization from metadata + if(metaDataJson.isMember("AnatomicRegionSequence")){ + frameAnatomyFG.getAnatomy().getAnatomicRegion().set( + metaDataJson["AnatomicRegionSequence"]["CodeValue"].asCString(), + metaDataJson["AnatomicRegionSequence"]["CodingSchemeDesignator"].asCString(), + metaDataJson["AnatomicRegionSequence"]["CodeMeaning"].asCString()); + } else { + frameAnatomyFG.getAnatomy().getAnatomicRegion().set("T-D0050", "SRT", "Tissue"); + } + + CHECK_COND(parametricMap->addForAllFrames(frameAnatomyFG)); + + return EXIT_SUCCESS; +} + +int ParametricMapObject::initializeRWVMFG() { + FGRealWorldValueMapping::RWVMItem* realWorldValueMappingItem = + new FGRealWorldValueMapping::RWVMItem(); + + realWorldValueMappingItem->setRealWorldValueSlope(metaDataJson["RealWorldValueSlope"].asFloat()); + realWorldValueMappingItem->setRealWorldValueIntercept(0); + + // Calculate intensity range - required + MinMaxCalculatorType::Pointer calculator = MinMaxCalculatorType::New(); + calculator->SetImage(itkImage); + calculator->Compute(); + + realWorldValueMappingItem->setRealWorldValueFirstValueMappedSigned(calculator->GetMinimum()); + realWorldValueMappingItem->setRealWorldValueLastValueMappedSigned(calculator->GetMaximum()); + + if(metaDataJson.isMember("MeasurementUnitsCode")){ + CodeSequenceMacro& unitsCodeDcmtk = realWorldValueMappingItem->getMeasurementUnitsCode(); + unitsCodeDcmtk.set(metaDataJson["MeasurementUnitsCode"]["CodeValue"].asCString(), + metaDataJson["MeasurementUnitsCode"]["CodingSchemeDesignator"].asCString(), + metaDataJson["MeasurementUnitsCode"]["CodeMeaning"].asCString()); + cout << "Measurements units initialized to " << + dcmqi::Helper::codeSequenceMacroToString(unitsCodeDcmtk); + + realWorldValueMappingItem->setLUTExplanation(metaDataJson["MeasurementUnitsCode"]["CodeMeaning"].asCString()); + realWorldValueMappingItem->setLUTLabel(metaDataJson["MeasurementUnitsCode"]["CodeValue"].asCString()); + } + + /* + if(metaDataJson.isMember("QuantityValueCode")){ + ContentItemMacro* item = initializeContentItemMacro(CodeSequenceMacro("G-C1C6", "SRT", "Quantity"), + dcmqi::Helper::jsonToCodeSequenceMacro(metaDataJson["QuantityValueCode"])); + realWorldValueMappingItem->getEntireQuantityDefinitionSequence().push_back(item); + }*/ + + ContentItemMacro* quantity = new ContentItemMacro; + CodeSequenceMacro* qCodeName = new CodeSequenceMacro("G-C1C6", "SRT", "Quantity"); + CodeSequenceMacro* qSpec = createCodeSequenceFromMetadata("QuantityValueCode"); + + if (!quantity || !qSpec || !qCodeName) + { + return EXIT_FAILURE; + } + + quantity->getEntireConceptNameCodeSequence().push_back(qCodeName); + quantity->getEntireConceptCodeSequence().push_back(qSpec); + realWorldValueMappingItem->getEntireQuantityDefinitionSequence().push_back(quantity); + quantity->setValueType(ContentItemMacro::VT_CODE); + + /* + // TODO: factor out defined CodeSequenceMacros into definitions as in dcmsr/include/dcmtk/dcmsr/codes + if(metaDataJson.isMember("MeasurementMethodCode")){ + ContentItemMacro* item = initializeContentItemMacro(CodeSequenceMacro("G-C306", "SRT", "Measurement Method"), + dcmqi::Helper::jsonToCodeSequenceMacro(metaDataJson["MeasurementMethodCode"])); + realWorldValueMappingItem->getEntireQuantityDefinitionSequence().push_back(item); + } */ + + if(metaDataJson.isMember("MeasurementMethodCode")){ + ContentItemMacro* measureMethod = new ContentItemMacro; + CodeSequenceMacro* qCodeName = new CodeSequenceMacro("G-C306", "SRT", "Measurement Method"); + CodeSequenceMacro* qSpec = createCodeSequenceFromMetadata("MeasurementMethodCode"); + + if (!measureMethod || !qSpec || !qCodeName) + { + return EXIT_FAILURE; + } + + measureMethod->getEntireConceptNameCodeSequence().push_back(qCodeName); + measureMethod->getEntireConceptCodeSequence().push_back(qSpec); + realWorldValueMappingItem->getEntireQuantityDefinitionSequence().push_back(measureMethod); + measureMethod->setValueType(ContentItemMacro::VT_CODE); + } + + /* + if(metaDataJson.isMember("ModelFittingMethodCode")){ + // TODO: update this once CP-1665 is integrated into the standard + ContentItemMacro* item = initializeContentItemMacro(CodeSequenceMacro("113241", "DCM", "Model Fitting Method"), + dcmqi::Helper::jsonToCodeSequenceMacro(metaDataJson["ModelFittingMethodCode"])); + realWorldValueMappingItem->getEntireQuantityDefinitionSequence().push_back(item); + } */ + + if(metaDataJson.isMember("ModelFittingMethodCode")){ + ContentItemMacro* fittingMethod = new ContentItemMacro; + CodeSequenceMacro* qCodeName = new CodeSequenceMacro("113241", "DCM", "Model fitting method"); + CodeSequenceMacro* qSpec = createCodeSequenceFromMetadata("ModelFittingMethodCode"); + + if (!fittingMethod || !qSpec || !qCodeName) + { + return EXIT_FAILURE; + } + + fittingMethod->getEntireConceptNameCodeSequence().push_back(qCodeName); + fittingMethod->getEntireConceptCodeSequence().push_back(qSpec); + realWorldValueMappingItem->getEntireQuantityDefinitionSequence().push_back(fittingMethod); + fittingMethod->setValueType(ContentItemMacro::VT_CODE); + } + + + if(metaDataJson.isMember("SourceImageDiffusionBValues")) { + for (int bvalId = 0; bvalId < metaDataJson["SourceImageDiffusionBValues"].size(); bvalId++) { + ContentItemMacro *bval = new ContentItemMacro; + CodeSequenceMacro *bvalUnits = new CodeSequenceMacro("s/mm2", "UCUM", "s/mm2"); + // TODO: update this once CP-1665 is integrated into the standard + CodeSequenceMacro *qCodeName = new CodeSequenceMacro("113240", "DCM", "Source image diffusion b-value"); + + if (!bval || !bvalUnits || !qCodeName) { + return EXIT_FAILURE; + } + + bval->setValueType(ContentItemMacro::VT_NUMERIC); + bval->getEntireConceptNameCodeSequence().push_back(qCodeName); + bval->getEntireMeasurementUnitsCodeSequence().push_back(bvalUnits); + if (bval->setNumericValue(metaDataJson["SourceImageDiffusionBValues"][bvalId].asCString()).bad()) + cout << "Failed to insert the value!" << endl;; + realWorldValueMappingItem->getEntireQuantityDefinitionSequence().push_back(bval); + } + } + + rwvmFG.getRealWorldValueMapping().push_back(realWorldValueMappingItem); + parametricMap->addForAllFrames(rwvmFG); + + return EXIT_SUCCESS; +} + +CodeSequenceMacro *ParametricMapObject::createCodeSequenceFromMetadata(const string &codeName) const { + return new CodeSequenceMacro( + metaDataJson[codeName]["CodeValue"].asCString(), + metaDataJson[codeName]["CodingSchemeDesignator"].asCString(), + metaDataJson[codeName]["CodeMeaning"].asCString()); +} + + +int ParametricMapObject::initializeFromDICOM(DcmDataset * sourceDataset) { + + sourceRepresentationType = DICOM_REPR; + dcmRepresentation = sourceDataset; + + DcmRLEDecoderRegistration::registerCodecs(); + + OFLogger dcemfinfLogger = OFLog::getLogger("qiicr.apps"); + dcemfinfLogger.setLogLevel(dcmtk::log4cplus::OFF_LOG_LEVEL); + + OFvariant result = DPMParametricMapIOD::loadDataset(*sourceDataset); + if (OFCondition* pCondition = OFget(&result)) { + throw -1; + } + + parametricMap = *OFget(&result); + + initializeMetaDataFromDICOM(); + createITKParametricMap(); + + return EXIT_SUCCESS; +} + +void ParametricMapObject::initializeMetaDataFromDICOM() { + + OFString temp; + parametricMap->getSeries().getSeriesDescription(temp); + metaDataJson["SeriesDescription"] = temp.c_str(); + + parametricMap->getSeries().getSeriesNumber(temp); + metaDataJson["SeriesNumber"] = temp.c_str(); + + parametricMap->getIODGeneralImageModule().getInstanceNumber(temp); + metaDataJson["InstanceNumber"] = temp.c_str(); + + using namespace dcmqi; + + parametricMap->getSeries().getBodyPartExamined(temp); + metaDataJson["BodyPartExamined"] = temp.c_str(); + + if (parametricMap->getNumberOfFrames() > 0) { + FGInterface& fg = parametricMap->getFunctionalGroups(); + FGRealWorldValueMapping* rw = OFstatic_cast(FGRealWorldValueMapping*, + fg.get(0, DcmFGTypes::EFG_REALWORLDVALUEMAPPING)); + if (rw->getRealWorldValueMapping().size() > 0) { + FGRealWorldValueMapping::RWVMItem *item = rw->getRealWorldValueMapping()[0]; + metaDataJson["MeasurementUnitsCode"] = Helper::codeSequence2Json(item->getMeasurementUnitsCode()); + + Float64 slope; + item->getData().findAndGetFloat64(DCM_RealWorldValueSlope, slope); + metaDataJson["RealWorldValueSlope"] = slope; + + vector diffusionBValues; + + for(int quantIdx=0; quantIdxgetEntireQuantityDefinitionSequence().size(); quantIdx++) { +// TODO: what if there are more than one? + ContentItemMacro* macro = item->getEntireQuantityDefinitionSequence()[quantIdx]; + CodeSequenceMacro* codeSequence= macro->getConceptNameCodeSequence(); + if (codeSequence != NULL) { + OFString codeMeaning; + codeSequence->getCodeMeaning(codeMeaning); + OFString designator, meaning, value; + + if (codeMeaning == "Quantity") { + CodeSequenceMacro* quantityValueCode = macro->getConceptCodeSequence(); + if (quantityValueCode != NULL) { + metaDataJson["QuantityValueCode"] = Helper::codeSequence2Json(*quantityValueCode); + } + } else if (codeMeaning == "Measurement Method") { + CodeSequenceMacro* measurementMethodValueCode = macro->getConceptCodeSequence(); + if (measurementMethodValueCode != NULL) { + metaDataJson["MeasurementMethodCode"] = Helper::codeSequence2Json(*measurementMethodValueCode); + } + } else if (codeMeaning == "Source image diffusion b-value") { + macro->getNumericValue(value); + diffusionBValues.push_back(value.c_str()); + } + } + } + + if (diffusionBValues.size() > 0) { + metaDataJson["SourceImageDiffusionBValues"] = Json::Value(Json::arrayValue); + for (vector::iterator it = diffusionBValues.begin() ; it != diffusionBValues.end(); ++it) + metaDataJson["SourceImageDiffusionBValues"].append(*it); + } + } + + FGDerivationImage* derivationImage = OFstatic_cast(FGDerivationImage*, fg.get(0, DcmFGTypes::EFG_DERIVATIONIMAGE)); + if(derivationImage) { + OFVector &derivationImageItems = derivationImage->getDerivationImageItems(); + if (derivationImageItems.size() > 0) { + DerivationImageItem *derivationImageItem = derivationImageItems[0]; + CodeSequenceMacro *derivationCode = derivationImageItem->getDerivationCodeItems()[0]; + if (derivationCode != NULL) { + metaDataJson["DerivationCode"] = Helper::codeSequence2Json(*derivationCode); + } + } + } + + FGFrameAnatomy* fa = OFstatic_cast(FGFrameAnatomy*, fg.get(0, DcmFGTypes::EFG_FRAMEANATOMY)); + metaDataJson["AnatomicRegionSequence"] = Helper::codeSequence2Json(fa->getAnatomy().getAnatomicRegion()); + + FGFrameAnatomy::LATERALITY frameLaterality; + fa->getLaterality(frameLaterality); + metaDataJson["FrameLaterality"] = fa->laterality2Str(frameLaterality).c_str(); + } +} + +bool ParametricMapObject::isDerivationFGRequired(vector >& slice2frame) { + // if there is a derivation item for at least one frame, DerivationImageSequence must be present for every frame. + unsigned nSlices = itkImage->GetLargestPossibleRegion().GetSize()[2]; + for (unsigned long sliceNumber=0;sliceNumber >& slice2frame){ + + FGPlanePosPatient* fgppp = FGPlanePosPatient::createMinimal("1","1","1"); + FGFrameContent* fgfc = new FGFrameContent(); + FGDerivationImage* fgder = new FGDerivationImage(); + OFVector perFrameFGs; + + unsigned nSlices = itkImage->GetLargestPossibleRegion().GetSize()[2]; + + bool derivationFGRequired = isDerivationFGRequired(slice2frame); + for (unsigned long sliceNumber = 0;sliceNumber < nSlices; sliceNumber++) { + + perFrameFGs.push_back(fgppp); + perFrameFGs.push_back(fgfc); + + fgder->clearData(); + if(!slice2frame[sliceNumber].empty()){ + if(metaDataJson.isMember("DerivationCode")){ + CodeSequenceMacro purposeOfReference = CodeSequenceMacro("121322","DCM","Source image for image processing operation"); + CodeSequenceMacro derivationCode = CodeSequenceMacro(metaDataJson["DerivationCode"]["CodeValue"].asCString(), + metaDataJson["DerivationCode"]["CodingSchemeDesignator"].asCString(), + metaDataJson["DerivationCode"]["CodeMeaning"].asCString()); + addDerivationItemToDerivationFG(fgder, slice2frame[sliceNumber], purposeOfReference, derivationCode); + } else { + addDerivationItemToDerivationFG(fgder, slice2frame[sliceNumber]); + } + } + + if(derivationFGRequired) + perFrameFGs.push_back(fgder); + + Float32ITKImageType::RegionType sliceRegion; + Float32ITKImageType::IndexType sliceIndex; + Float32ITKImageType::SizeType inputSize = itkImage->GetBufferedRegion().GetSize(); + + sliceIndex[0] = 0; + sliceIndex[1] = 0; + sliceIndex[2] = sliceNumber; + + inputSize[2] = 1; + + sliceRegion.SetIndex(sliceIndex); + sliceRegion.SetSize(inputSize); + + const unsigned frameSize = inputSize[0] * inputSize[1]; + + OFVector data(frameSize); + + itk::ImageRegionConstIteratorWithIndex sliceIterator(itkImage, sliceRegion); + + unsigned framePixelCnt = 0; + for(sliceIterator.GoToBegin();!sliceIterator.IsAtEnd(); ++sliceIterator, ++framePixelCnt){ + data[framePixelCnt] = sliceIterator.Get(); + Float32ITKImageType::IndexType idx = sliceIterator.GetIndex(); + // cout << framePixelCnt << " " << idx[1] << "," << idx[0] << endl; + } + + // Plane Position + Float32ITKImageType::PointType sliceOriginPoint; + itkImage->TransformIndexToPhysicalPoint(sliceIndex, sliceOriginPoint); + fgppp->setImagePositionPatient( + dcmqi::Helper::floatToStrScientific(sliceOriginPoint[0]).c_str(), + dcmqi::Helper::floatToStrScientific(sliceOriginPoint[1]).c_str(), + dcmqi::Helper::floatToStrScientific(sliceOriginPoint[2]).c_str()); + + // Frame Content + OFCondition result = fgfc->setDimensionIndexValues(sliceNumber+1 /* value within dimension */, 0 /* first dimension */); + + DPMParametricMapIOD::FramesType frames = parametricMap->getFrames(); + result = OFget >(&frames)->addFrame(&*data.begin(), frameSize, perFrameFGs); + + perFrameFGs.clear(); + + } + + return EXIT_SUCCESS; +} diff --git a/libsrc/SegmentVolume.cpp b/libsrc/SegmentVolume.cpp new file mode 100644 index 00000000..e37c473a --- /dev/null +++ b/libsrc/SegmentVolume.cpp @@ -0,0 +1,5 @@ +// +// Created by Andrey Fedorov on 3/9/17. +// + +#include "dcmqi/SegmentVolume.h" diff --git a/libsrc/ImageSEGConverter.cpp b/libsrc/SegmentationImageConverter.cpp similarity index 95% rename from libsrc/ImageSEGConverter.cpp rename to libsrc/SegmentationImageConverter.cpp index 1eaf4591..a294d3fb 100644 --- a/libsrc/ImageSEGConverter.cpp +++ b/libsrc/SegmentationImageConverter.cpp @@ -1,11 +1,40 @@ // DCMQI includes -#include "dcmqi/ImageSEGConverter.h" +#include "dcmqi/SegmentationImageConverter.h" +#include "dcmqi/SegmentationImageObject.h" +using namespace std; namespace dcmqi { - DcmDataset* ImageSEGConverter::itkimage2dcmSegmentation(vector dcmDatasets, + DcmDataset* itkimage2paramapReplacement(vector dcmDatasets, + vector segmentations, + const string &metaData, + bool skipEmptySlices) { + SegmentationImageObject seg; + seg.initializeFromITK(segmentations, metaData, dcmDatasets, skipEmptySlices); + + DcmDataset *output = new DcmDataset(); + seg.getDICOMRepresentation(*output); + + return output; + } + + pair , string> dcmSegmentation2itkimageReplacement(DcmDataset *segDataset) { + + SegmentationImageObject seg; + seg.initializeFromDICOM(segDataset); + + Json::StyledWriter styledWriter; + std::stringstream ss; + + ss << styledWriter.write(seg.getMetaDataJson()); + + return pair , string>(seg.getITKRepresentation(), + ss.str()); + }; + + DcmDataset* SegmentationImageConverter::itkimage2dcmSegmentation(vector dcmDatasets, vector segmentations, const string &metaData, bool skipEmptySlices) { @@ -116,7 +145,6 @@ namespace dcmqi { if((*vI).size()>0) hasDerivationImages = true; - FGPlanePosPatient* fgppp = FGPlanePosPatient::createMinimal("1","1","1"); FGFrameContent* fgfc = new FGFrameContent(); FGDerivationImage* fgder = new FGDerivationImage(); @@ -448,7 +476,7 @@ namespace dcmqi { } - pair , string> ImageSEGConverter::dcmSegmentation2itkimage(DcmDataset *segDataset) { + pair , string> SegmentationImageConverter::dcmSegmentation2itkimage(DcmDataset *segDataset) { DcmRLEDecoderRegistration::registerCodecs(); @@ -704,7 +732,7 @@ namespace dcmqi { return pair , string>(segment2image, metaInfo.getJSONOutputAsString()); } - void ImageSEGConverter::populateMetaInformationFromDICOM(DcmDataset *segDataset, DcmSegmentation *segdoc, + void SegmentationImageConverter::populateMetaInformationFromDICOM(DcmDataset *segDataset, DcmSegmentation *segdoc, JSONSegmentationMetaInformationHandler &metaInfo) { OFString creatorName, sessionID, timePointID, seriesDescription, seriesNumber, instanceNumber, bodyPartExamined, coordinatingCenter; diff --git a/libsrc/SegmentationImageObject.cpp b/libsrc/SegmentationImageObject.cpp new file mode 100644 index 00000000..5b7d6be0 --- /dev/null +++ b/libsrc/SegmentationImageObject.cpp @@ -0,0 +1,714 @@ +// +// Created by Andrey Fedorov on 3/11/17. +// + +#include "dcmqi/SegmentationImageObject.h" + +int SegmentationImageObject::initializeFromITK(vector itkSegmentations, + const string &metaDataStr, vector derivationDatasets, + bool skipEmptySlices) { + + setDerivationDatasets(derivationDatasets); + + sourceRepresentationType = ITK_REPR; + + // TODO: think about the following code. Do all input segmentations have the same extent? + itkImages = itkSegmentations; + + initializeMetaDataFromString(metaDataStr); + +// if(!metaDataIsComplete()){ +// updateMetaDataFromDICOM(derivationDatasets); +// } + + initializeVolumeGeometry(); + + // NB: the sequence of steps initializing different components of the object parallels that + // in the original converter function. It probably makes sense to revisit the sequence + // of these steps. It does not necessarily need to happen in this order. + initializeEquipmentInfo(); + initializeContentIdentification(); + +// OFString frameOfRefUID; +// if(!segdoc->getFrameOfReference().getFrameOfReferenceUID(frameOfRefUID).good()){ +// // TODO: add FoR UID to the metadata JSON and check that before generating one! +// char frameOfRefUIDchar[128]; +// dcmGenerateUniqueIdentifier(frameOfRefUIDchar, QIICR_UID_ROOT); +// CHECK_COND(segdoc->getFrameOfReference().setFrameOfReferenceUID(frameOfRefUIDchar)); +// } +// + + createDICOMSegmentation(); + + // populate metadata about patient/study, from derivation + // datasets or from metadata + initializeCompositeContext(); + + initializeSeriesSpecificAttributes(segmentation->getSeries(), + segmentation->getGeneralImage()); + + initializePixelMeasuresFG(); + CHECK_COND(segmentation->addForAllFrames(pixelMeasuresFG)); + + initializePlaneOrientationFG(); + CHECK_COND(segmentation->addForAllFrames(planeOrientationPatientFG)); + + vector > slice2frame; + mapVolumeSlicesToDICOMFrames(derivationDatasets, slice2frame); + + // NB: this assumes all segmentation files have the same dimensions; alternatively, need to + // do this operation for each segmentation file + vector > slice2derimg = getSliceMapForSegmentation2DerivationImage(derivationDatasets, + itkSegmentations[0]); + + initializeCommonInstanceReferenceModule(segmentation->getCommonInstanceReference(), slice2frame); + + // initial frames ... + + FGPlanePosPatient* fgppp = FGPlanePosPatient::createMinimal("1","1","1"); + FGFrameContent* fgfc = new FGFrameContent(); + FGDerivationImage* fgder = new FGDerivationImage(); + OFVector perFrameFGs; + + perFrameFGs.push_back(fgppp); + perFrameFGs.push_back(fgfc); + if(hasDerivationImages(slice2derimg)) + perFrameFGs.push_back(fgder); + + int uidfound = 0, uidnotfound = 0; + const unsigned frameSize = volumeGeometry.extent[0] * volumeGeometry.extent[1]; + Uint8 *frameData = new Uint8[frameSize]; + +// for(size_t segFileNumber=0; segFileNumberSetInput(itkImages[segFileNumber]); +// l2lm->Update(); +// +// typedef LabelToLabelMapFilterType::OutputImageType::LabelObjectType LabelType; +// typedef itk::LabelStatisticsImageFilter LabelStatisticsType; +// +// LabelStatisticsType::Pointer labelStats = LabelStatisticsType::New(); +// +// cout << "Found " << l2lm->GetOutput()->GetNumberOfLabelObjects() << " label(s)" << endl; +// labelStats->SetInput(itkImages[segFileNumber]); +// labelStats->SetLabelInput(itkImages[segFileNumber]); +// labelStats->Update(); +// +// // TODO: implement crop segments box +// +// for(unsigned segLabelNumber=0 ; segLabelNumberGetOutput()->GetNumberOfLabelObjects();segLabelNumber++){ +// LabelType* labelObject = l2lm->GetOutput()->GetNthLabelObject(segLabelNumber); +// short label = labelObject->GetLabel(); +// +// if(!label){ +// cout << "Skipping label 0" << endl; +// continue; +// } +// +// cout << "Processing label " << label << endl; +// +// LabelStatisticsType::BoundingBoxType bbox = labelStats->GetBoundingBox(label); +// unsigned firstSlice, lastSlice; +// //bool skipEmptySlices = true; // TODO: what to do with that line? +// //bool skipEmptySlices = false; // TODO: what to do with that line? +// if(skipEmptySlices){ +// firstSlice = bbox[4]; +// lastSlice = bbox[5]+1; +// } else { +// firstSlice = 0; +// lastSlice = inputSize[2]; +// } +// +// cout << "Total non-empty slices that will be encoded in SEG for label " << +// label << " is " << lastSlice-firstSlice << endl << +// " (inclusive from " << firstSlice << " to " << +// lastSlice << ")" << endl; +// +// DcmSegment* segment = NULL; +// if(metaInfo.segmentsAttributesMappingList[segFileNumber].find(label) == metaInfo.segmentsAttributesMappingList[segFileNumber].end()){ +// cerr << "ERROR: Failed to match label from image to the segment metadata!" << endl; +// return NULL; +// } +// +// SegmentAttributes* segmentAttributes = metaInfo.segmentsAttributesMappingList[segFileNumber][label]; +// +// DcmSegTypes::E_SegmentAlgoType algoType = DcmSegTypes::SAT_UNKNOWN; +// string algoName = ""; +// string algoTypeStr = segmentAttributes->getSegmentAlgorithmType(); +// if(algoTypeStr == "MANUAL"){ +// algoType = DcmSegTypes::SAT_MANUAL; +// } else { +// if(algoTypeStr == "AUTOMATIC") +// algoType = DcmSegTypes::SAT_AUTOMATIC; +// if(algoTypeStr == "SEMIAUTOMATIC") +// algoType = DcmSegTypes::SAT_SEMIAUTOMATIC; +// +// algoName = segmentAttributes->getSegmentAlgorithmName(); +// if(algoName == ""){ +// cerr << "ERROR: Algorithm name must be specified for non-manual algorithm types!" << endl; +// return NULL; +// } +// } +// +// CodeSequenceMacro* typeCode = segmentAttributes->getSegmentedPropertyTypeCodeSequence(); +// CodeSequenceMacro* categoryCode = segmentAttributes->getSegmentedPropertyCategoryCodeSequence(); +// assert(typeCode != NULL && categoryCode!= NULL); +// OFString segmentLabel; +// CHECK_COND(typeCode->getCodeMeaning(segmentLabel)); +// CHECK_COND(DcmSegment::create(segment, segmentLabel, *categoryCode, *typeCode, algoType, algoName.c_str())); +// +// if(segmentAttributes->getSegmentDescription().length() > 0) +// segment->setSegmentDescription(segmentAttributes->getSegmentDescription().c_str()); +// +// CodeSequenceMacro* typeModifierCode = segmentAttributes->getSegmentedPropertyTypeModifierCodeSequence(); +// if (typeModifierCode != NULL) { +// OFVector& modifiersVector = segment->getSegmentedPropertyTypeModifierCode(); +// modifiersVector.push_back(typeModifierCode); +// } +// +// GeneralAnatomyMacro &anatomyMacro = segment->getGeneralAnatomyCode(); +// if (segmentAttributes->getAnatomicRegionSequence() != NULL){ +// OFVector& anatomyMacroModifiersVector = anatomyMacro.getAnatomicRegionModifier(); +// CodeSequenceMacro& anatomicRegionSequence = anatomyMacro.getAnatomicRegion(); +// anatomicRegionSequence = *segmentAttributes->getAnatomicRegionSequence(); +// +// if(segmentAttributes->getAnatomicRegionModifierSequence() != NULL){ +// CodeSequenceMacro* anatomicRegionModifierSequence = segmentAttributes->getAnatomicRegionModifierSequence(); +// anatomyMacroModifiersVector.push_back(anatomicRegionModifierSequence); +// } +// } +// +// unsigned* rgb = segmentAttributes->getRecommendedDisplayRGBValue(); +// unsigned cielabScaled[3]; +// float cielab[3], ciexyz[3]; +// +// Helper::getCIEXYZFromRGB(&rgb[0],&ciexyz[0]); +// Helper::getCIELabFromCIEXYZ(&ciexyz[0],&cielab[0]); +// Helper::getIntegerScaledCIELabFromCIELab(&cielab[0],&cielabScaled[0]); +// CHECK_COND(segment->setRecommendedDisplayCIELabValue(cielabScaled[0],cielabScaled[1],cielabScaled[2])); +// +// Uint16 segmentNumber; +// CHECK_COND(segdoc->addSegment(segment, segmentNumber /* returns logical segment number */)); +// +// // TODO: make it possible to skip empty frames (optional) +// // iterate over slices for an individual label and populate output frames +// for(unsigned sliceNumber=firstSlice;sliceNumbersetStackID("1"); // all frames go into the same stack +// CHECK_COND(fgfc->setDimensionIndexValues(segmentNumber, 0)); +// CHECK_COND(fgfc->setDimensionIndexValues(sliceNumber-firstSlice+1, 1)); +// //ostringstream inStackPosSStream; // StackID is not present/needed +// //inStackPosSStream << s+1; +// //fracon->setInStackPositionNumber(s+1); +// +// // PerFrame FG: PlanePositionSequence +// { +// ShortImageType::PointType sliceOriginPoint; +// ShortImageType::IndexType sliceOriginIndex; +// sliceOriginIndex.Fill(0); +// sliceOriginIndex[2] = sliceNumber; +// itkImages[segFileNumber]->TransformIndexToPhysicalPoint(sliceOriginIndex, sliceOriginPoint); +// ostringstream pppSStream; +// if(sliceNumber>0){ +// ShortImageType::PointType prevOrigin; +// ShortImageType::IndexType prevIndex; +// prevIndex.Fill(0); +// prevIndex[2] = sliceNumber-1; +// itkImages[segFileNumber]->TransformIndexToPhysicalPoint(prevIndex, prevOrigin); +// } +// fgppp->setImagePositionPatient( +// Helper::floatToStrScientific(sliceOriginPoint[0]).c_str(), +// Helper::floatToStrScientific(sliceOriginPoint[1]).c_str(), +// Helper::floatToStrScientific(sliceOriginPoint[2]).c_str()); +// } +// +// /* Add frame that references this segment */ +// { +// ShortImageType::RegionType sliceRegion; +// ShortImageType::IndexType sliceIndex; +// ShortImageType::SizeType sliceSize; +// +// sliceIndex[0] = 0; +// sliceIndex[1] = 0; +// sliceIndex[2] = sliceNumber; +// +// sliceSize[0] = volumeGeometry.extent[0]; +// sliceSize[1] = volumeGeometry.extent[1]; +// sliceSize[2] = 1; +// +// sliceRegion.SetIndex(sliceIndex); +// sliceRegion.SetSize(sliceSize); +// +// unsigned framePixelCnt = 0; +// itk::ImageRegionConstIteratorWithIndex sliceIterator(itkImages[segFileNumber], sliceRegion); +// for(sliceIterator.GoToBegin();!sliceIterator.IsAtEnd();++sliceIterator,++framePixelCnt){ +// if(sliceIterator.Get() == label){ +// frameData[framePixelCnt] = 1; +// ShortImageType::IndexType idx = sliceIterator.GetIndex(); +// //cout << framePixelCnt << " " << idx[1] << "," << idx[0] << endl; +// } else +// frameData[framePixelCnt] = 0; +// } +// +// /* +// if(sliceNumber>=dcmDatasets.size()){ +// cerr << "ERROR: trying to access missing DICOM Slice! And sorry, multi-frame not supported at the moment..." << endl; +// return NULL; +// }*/ +// +// OFVector siVector; +// for(size_t derImageInstanceNum=0; +// derImageInstanceNum0){ +// +// DerivationImageItem *derimgItem; +// CHECK_COND(fgder->addDerivationImageItem(CodeSequenceMacro("113076","DCM","Segmentation"),"",derimgItem)); +// +// cout << "Total of " << siVector.size() << " source image items will be added" << endl; +// +// OFVector srcimgItems; +// CHECK_COND(derimgItem->addSourceImageItems(siVector, +// CodeSequenceMacro("121322","DCM","Source image for image processing operation"), +// srcimgItems)); +// +// if(1){ +// // initialize class UID and series instance UID +// ImageSOPInstanceReferenceMacro &instRef = srcimgItems[0]->getImageSOPInstanceReference(); +// OFString instanceUID; +// CHECK_COND(instRef.getReferencedSOPClassUID(classUID)); +// CHECK_COND(instRef.getReferencedSOPInstanceUID(instanceUID)); +// +// if(instanceUIDs.find(instanceUID) == instanceUIDs.end()){ +// SOPInstanceReferenceMacro *refinstancesItem = new SOPInstanceReferenceMacro(); +// CHECK_COND(refinstancesItem->setReferencedSOPClassUID(classUID)); +// CHECK_COND(refinstancesItem->setReferencedSOPInstanceUID(instanceUID)); +// refinstances.push_back(refinstancesItem); +// instanceUIDs.insert(instanceUID); +// uidnotfound++; +// } else { +// uidfound++; +// } +// } +// } +// +// CHECK_COND(segdoc->addFrame(frameData, segmentNumber, perFrameFGs)); +// +// // remove derivation image FG from the per-frame FGs, only if applicable! +// if(siVector.size()>0){ +// // clean up for the next frame +// fgder->clearData(); +// } +// +// } +// } +// } +// } +// +// delete fgfc; +// delete fgppp; +// delete fgder; + + return EXIT_SUCCESS; + +} + +bool SegmentationImageObject::hasDerivationImages(vector > &slice2derimg) const { + for (vector >::const_iterator vI = slice2derimg.begin(); vI != slice2derimg.end(); ++vI) { + if ((*vI).size() > 0) { + return true; + } + } + return false; +} + +int SegmentationImageObject::initializeEquipmentInfo() { + if(sourceRepresentationType == ITK_REPR){ + generalEquipmentInfoModule = IODGeneralEquipmentModule::EquipmentInfo(QIICR_MANUFACTURER, + QIICR_DEVICE_SERIAL_NUMBER, + QIICR_MANUFACTURER_MODEL_NAME, + QIICR_SOFTWARE_VERSIONS); + } else { // DICOM_REPR + } + return EXIT_SUCCESS; +} + +int SegmentationImageObject::initializeVolumeGeometry() { + if(sourceRepresentationType == ITK_REPR){ + + ShortToDummyCasterType::Pointer caster = + ShortToDummyCasterType::New(); + // right now assuming that all input segmentations have the same volume extent + // TODO: make this work when input segmentation have different geometry + caster->SetInput(itkImages[0]); + caster->Update(); + + MultiframeObject::initializeVolumeGeometryFromITK(caster->GetOutput()); + } else { + MultiframeObject::initializeVolumeGeometryFromDICOM(segmentation->getFunctionalGroups()); + } + return EXIT_SUCCESS; +} + +int SegmentationImageObject::createDICOMSegmentation() { + + DcmSegmentation::createBinarySegmentation(segmentation, + volumeGeometry.extent[1], + volumeGeometry.extent[0], + generalEquipmentInfoModule, + contentIdentificationMacro); // content identification + + /* Initialize dimension module */ + IODMultiframeDimensionModule &mfdim = segmentation->getDimensions(); + OFString dimUID = dcmqi::Helper::generateUID(); + CHECK_COND(mfdim.addDimensionIndex(DCM_ReferencedSegmentNumber, dimUID, DCM_SegmentIdentificationSequence, + DcmTag(DCM_ReferencedSegmentNumber).getTagName())); + CHECK_COND(mfdim.addDimensionIndex(DCM_ImagePositionPatient, dimUID, + DCM_PlanePositionSequence, "ImagePositionPatient")); + + return EXIT_SUCCESS; +} + +int SegmentationImageObject::initializeCompositeContext() { + if(derivationDcmDatasets.size()){ + CHECK_COND(segmentation->import(*derivationDcmDatasets[0], + OFTrue, // Patient + OFTrue, // Study + OFTrue, // Frame of reference + OFFalse)); // Series + } else { + // TODO: once we support passing of composite context in metadata, propagate it into segmentation here + return EXIT_FAILURE; + } + return EXIT_SUCCESS; +} + +int SegmentationImageObject::initializeFromDICOM(DcmDataset* sourceDataset) { + + dcmRepresentation = sourceDataset; + +// TODO: add SegmentationImageObject to namespace dcmqi + using namespace dcmqi; + + DcmRLEDecoderRegistration::registerCodecs(); + + OFLogger dcemfinfLogger = OFLog::getLogger("qiicr.apps"); + dcemfinfLogger.setLogLevel(dcmtk::log4cplus::OFF_LOG_LEVEL); + + OFCondition cond = DcmSegmentation::loadDataset(*sourceDataset, segmentation); + if(!segmentation){ + cerr << "Failed to load seg! " << cond.text() << endl; + throw -1; + } + + initializeVolumeGeometryFromDICOM(segmentation->getFunctionalGroups()); + itkImage = volumeGeometry.getITKRepresentation(); + vector< pair > matchedFramesWithSlices = matchFramesWithSegmentIdAndSliceNumber( + segmentation->getFunctionalGroups()); + unpackFramesAndWriteSegmentImage(matchedFramesWithSlices); + initializeMetaDataFromDICOM(sourceDataset); + + return EXIT_SUCCESS; +} + +vector< pair > SegmentationImageObject::matchFramesWithSegmentIdAndSliceNumber(FGInterface &fgInterface) { + + vector< pair > matchedFramesWithSlices; + + for(size_t frameId=0;frameIdgetImagePositionPatient(planposStr, j).good()){ + frameOriginPoint[j] = atof(planposStr.c_str()); + } + } + clog << "Image size: " << segment2image[segmentId]->GetBufferedRegion().GetSize() << endl; + + if(!segment2image[segmentId]->TransformPhysicalPointToIndex(frameOriginPoint, frameOriginIndex)){ + cerr << "ERROR: Frame " << frameId << " origin " << frameOriginPoint << + " is outside image geometry!" << frameOriginIndex << endl; + cerr << "Image size: " << segment2image[segmentId]->GetBufferedRegion().GetSize() << endl; + throw -1; + } + matchedFramesWithSlices.push_back(make_pair(segmentId, frameOriginIndex[2])); + } + return matchedFramesWithSlices; +} + +int SegmentationImageObject::unpackFramesAndWriteSegmentImage( + vector< pair > matchingSegmentIDsAndSliceNumbers) { + + for (std::vector< pair >::iterator it = matchingSegmentIDsAndSliceNumbers.begin(); + it != matchingSegmentIDsAndSliceNumbers.end(); ++it) { + + unsigned frameId = it - matchingSegmentIDsAndSliceNumbers.begin(); + Uint16 segmentId = (*it).first; + long sliceNumber = (*it).second; + + const DcmIODTypes::Frame *frame = segmentation->getFrame(frameId); + + DcmIODTypes::Frame *unpackedFrame = NULL; + + if (segmentation->getSegmentationType() == DcmSegTypes::ST_BINARY) { + unpackedFrame = DcmSegUtils::unpackBinaryFrame(frame, + volumeGeometry.extent[1], // Rows + volumeGeometry.extent[0]); // Cols + } else { + unpackedFrame = new DcmIODTypes::Frame(*frame); + } + + for (unsigned row = 0; row < volumeGeometry.extent[1]; row++) { + for (unsigned col = 0; col < volumeGeometry.extent[0]; col++) { + ShortImageType::PixelType pixel; + unsigned bitCnt = row * volumeGeometry.extent[0] + col; + pixel = unpackedFrame->pixData[bitCnt]; + + if (pixel != 0) { + ShortImageType::IndexType index; + index[0] = col; + index[1] = row; + index[2] = sliceNumber; + segment2image[segmentId]->SetPixel(index, segmentId); + } + } + } + + if (unpackedFrame != NULL) + delete unpackedFrame; + } + + return EXIT_SUCCESS; +} + +int SegmentationImageObject::createNewSegmentImage(Uint16 segmentId) { + typedef itk::ImageDuplicator DuplicatorType; + DuplicatorType::Pointer dup = DuplicatorType::New(); + dup->SetInputImage(itkImage); + dup->Update(); + ShortImageType::Pointer newSegmentImage = dup->GetOutput(); + newSegmentImage->FillBuffer(0); + segment2image[segmentId] = newSegmentImage; + return EXIT_SUCCESS; +} + +int SegmentationImageObject::initializeMetaDataFromDICOM(DcmDataset *segDataset) { + + OFString temp; + segmentation->getSeries().getSeriesDescription(temp); + metaDataJson["SeriesDescription"] = temp.c_str(); + + segmentation->getSeries().getSeriesNumber(temp); + metaDataJson["SeriesNumber"] = temp.c_str(); + + segDataset->findAndGetOFString(DCM_InstanceNumber, temp); + metaDataJson["InstanceNumber"] = temp.c_str(); + + segmentation->getSeries().getBodyPartExamined(temp); + metaDataJson["BodyPartExamined"] = temp.c_str(); + + segmentation->getContentIdentification().getContentCreatorName(temp); + metaDataJson["ContentCreatorName"] = temp.c_str(); + + segDataset->findAndGetOFString(DCM_ClinicalTrialTimePointID, temp); + metaDataJson["ClinicalTrialTimePointID"] = temp.c_str(); + + segDataset->findAndGetOFString(DCM_ClinicalTrialSeriesID, temp); + metaDataJson["ClinicalTrialSeriesID"] = temp.c_str(); + + segDataset->findAndGetOFString(DCM_ClinicalTrialCoordinatingCenterName, temp); + if (temp.size()) + metaDataJson["ClinicalTrialCoordinatingCenterName"] = temp.c_str(); + + metaDataJson["segmentAttributes"] = getSegmentAttributesMetadata(); + + return EXIT_SUCCESS; +} + +Json::Value SegmentationImageObject::getSegmentAttributesMetadata() { + + using namespace dcmqi; + + FGInterface &fgInterface = segmentation->getFunctionalGroups(); + + Json::Value values(Json::arrayValue); + vector processedSegmentIDs; + + for(size_t frameId=0;frameIdgetSegment(segmentId); + if (segment == NULL) { + cerr << "Failed to get segment for segment ID " << segmentId << endl; + continue; + } + + if(std::find(processedSegmentIDs.begin(), processedSegmentIDs.end(), segmentId) != processedSegmentIDs.end()) { + continue; + } + processedSegmentIDs.push_back(segmentId); + + // get CIELab color for the segment + Uint16 ciedcm[3]; + unsigned cielabScaled[3]; + float cielab[3], ciexyz[3]; + unsigned rgb[3]; + if (segment->getRecommendedDisplayCIELabValue( + ciedcm[0], ciedcm[1], ciedcm[2] + ).bad()) { + // NOTE: if the call above fails, it overwrites the values anyway, + // not sure if this is a dcmtk bug or not + ciedcm[0] = 43803; + ciedcm[1] = 26565; + ciedcm[2] = 37722; + cerr << "Failed to get CIELab values - initializing to default " << + ciedcm[0] << "," << ciedcm[1] << "," << ciedcm[2] << endl; + } + cielabScaled[0] = unsigned(ciedcm[0]); + cielabScaled[1] = unsigned(ciedcm[1]); + cielabScaled[2] = unsigned(ciedcm[2]); + + Helper::getCIELabFromIntegerScaledCIELab(&cielabScaled[0], &cielab[0]); + Helper::getCIEXYZFromCIELab(&cielab[0], &ciexyz[0]); + Helper::getRGBFromCIEXYZ(&ciexyz[0], &rgb[0]); + + Json::Value segmentEntry; + + OFString temp; + + segmentEntry["labelID"] = segmentId; + + segment->getSegmentDescription(temp); + segmentEntry["SegmentDescription"] = temp.c_str(); + + + DcmSegTypes::E_SegmentAlgoType algorithmType = segment->getSegmentAlgorithmType(); + string readableAlgorithmType = DcmSegTypes::algoType2OFString(algorithmType).c_str(); + segmentEntry["SegmentAlgorithmType"] = readableAlgorithmType; + + if (algorithmType == DcmSegTypes::SAT_UNKNOWN) { + cerr << "AlgorithmType is not valid with value " << readableAlgorithmType << endl; + throw -1; + } + if (algorithmType != DcmSegTypes::SAT_MANUAL) { + segment->getSegmentAlgorithmName(temp); + if (temp.length() > 0) + segmentEntry["SegmentAlgorithmName"] = temp.c_str(); + } + + Json::Value rgbArray(Json::arrayValue); + rgbArray.append(rgb[0]); + rgbArray.append(rgb[1]); + rgbArray.append(rgb[2]); + segmentEntry["recommendedDisplayRGBValue"] = rgbArray; + + segmentEntry["SegmentedPropertyCategoryCodeSequence"] = + Helper::codeSequence2Json(segment->getSegmentedPropertyCategoryCode()); + + segmentEntry["SegmentedPropertyTypeCodeSequence"] = + Helper::codeSequence2Json(segment->getSegmentedPropertyTypeCode()); + + if (segment->getSegmentedPropertyTypeModifierCode().size() > 0) { + segmentEntry["SegmentedPropertyTypeModifierCodeSequence"] = + Helper::codeSequence2Json(*(segment->getSegmentedPropertyTypeModifierCode())[0]); + } + + GeneralAnatomyMacro &anatomyMacro = segment->getGeneralAnatomyCode(); + CodeSequenceMacro &anatomicRegionSequence = anatomyMacro.getAnatomicRegion(); + if (anatomicRegionSequence.check(true).good()) { + segmentEntry["AnatomicRegionSequence"] = Helper::codeSequence2Json(anatomyMacro.getAnatomicRegion()); + } + if (anatomyMacro.getAnatomicRegionModifier().size() > 0) { + segmentEntry["AnatomicRegionModifierSequence"] = + Helper::codeSequence2Json(*(anatomyMacro.getAnatomicRegionModifier()[0])); + } + + Json::Value innerList(Json::arrayValue); + innerList.append(segmentEntry); + values.append(innerList); + } + return values; +} + +Uint16 SegmentationImageObject::getSegmentId(FGInterface &fgInterface, size_t frameId) const { + bool isPerFrame; + FGSegmentation *fgseg = + OFstatic_cast(FGSegmentation*,fgInterface.get(frameId, DcmFGTypes::EFG_SEGMENTATION, isPerFrame)); + assert(fgseg); + + Uint16 segmentId = -1; + if(fgseg->getReferencedSegmentNumber(segmentId).bad()){ + cerr << "Failed to get seg number!"; + throw -1; + } + return segmentId; +} + +template +vector > SegmentationImageObject::getSliceMapForSegmentation2DerivationImage(const vector dcmDatasets, + const ImageTypePointer &labelImage) { + // Find mapping from the segmentation slice number to the derivation image + // Assume that orientation of the segmentation is the same as the source series + unsigned numLabelSlices = labelImage->GetLargestPossibleRegion().GetSize()[2]; + vector > slice2derimg(numLabelSlices); + for(size_t i=0;ifindAndGetOFString(DCM_ImagePositionPatient, ippStr, j)); + ippPoint[j] = atof(ippStr.c_str()); + } + // NB: this will map slice origin to index without failure, unless the point is out + // of FOV bounds! + // TODO: do a better job matching volume slices by considering comparison of the origin + // and orientation of the slice within tolerance + if(!labelImage->TransformPhysicalPointToIndex(ippPoint, ippIndex)){ + //cout << "image position: " << ippPoint << endl; + //cerr << "ippIndex: " << ippIndex << endl; + // if certain DICOM instance does not map to a label slice, just skip it + continue; + } + slice2derimg[ippIndex[2]].push_back(i); + } + return slice2derimg; +} \ No newline at end of file