Snippets
Displaying: model_mri_dicom.cpp | View full page#include "stdafx.h" #include "CharLS/interface.h" #include <float.h> #include <algorithm> //section/table/figure/etc. references in this module apply to http://medical.nema.org/Dicom/2011/11_05pu.pdf //"Digital Imaging and Communications in Medicine (DICOM) // Part 5: Data Structures and Encoding" //...unless otherwise explicitly noted. (some things also explicitly refer to part 3) dicomOpts_t *g_dicomOpts = NULL; /* Private UID format (see Annex B): 1.2.840.xxxxx.3.152.235.2.12.187636473" root suffix In this example, the root is: 1 Identifies ISO 2 Identifies ANSI Member Body 840 Country code of a specific Member Body (U.S. for ANSI) xxxxx Identifies a specific Organization.(provided by ANSI) In this example the first two components of the suffix relate to the identification of the device: 3 Manufacturer or user defined device type 152 Manufacturer or user defined serial number The remaining four components of the suffix relate to the identification of the image: 235 Study number 2 Series number 12 Image number 187636473 Encoded date and time stamp of image acquisition */ enum EDicomTransferSyntax { kDTS_Unknown = 0, kDTS_ImplicitVRLittleEndian, kDTS_ExplicitVRLittleEndian, kDTS_ExplicitVRBigEndian, kDTS_RLE, //is_jpeg_transfer_syntax reflects kDTS_JpegBaseline as first jpeg transfer syntax kDTS_JpegBaseline, kDTS_JpegExtended, kDTS_JpegSpectral, kDTS_JpegProgressive, kDTS_JpegLossless, kDTS_JpegLosslessFirstOrder, kDTS_JpegLSLossless, kDTS_JpegLSNearLossless, kDTS_Jpeg2000LosslessReversible, kDTS_Jpeg2000LosslessOrLossy, kDTS_Jpeg2000Part2LosslessReversible, kDTS_Jpeg2000Part2LosslessOrLossy, //is_jpeg_transfer_syntax reflects kDTS_Jpeg2000Part2LosslessOrLossy as last jpeg transfer syntax kDTS_MPEG2_ML, kDTS_MPEG2_HL, kDTS_MPEG4, kDTS_MPEG4Temporal, kDTS_Deflate, kDTS_Count }; static const char *skTransferSyntaxStrings[kDTS_Count] = { "Unknown", //kDTS_Unknown "1.2.840.10008.1.2", //kDTS_ImplicitVRLittleEndian (A.1) "1.2.840.10008.1.2.1", //kDTS_ExplicitVRLittleEndian (A.2) "1.2.840.10008.1.2.2", //kDTS_ExplicitVRBigEndian (A.3) "1.2.840.10008.1.2.5", //kDTS_RLE (A.4.2) "1.2.840.10008.1.2.4.50", //kDTS_JpegBaseline (table A.4-3) "1.2.840.10008.1.2.4.51", //kDTS_JpegExtended (table A.4-3) "1.2.840.10008.1.2.4.53", //kDTS_JpegSpectral (not listed) "1.2.840.10008.1.2.4.55", //kDTS_JpegProgressive (not listed) "1.2.840.10008.1.2.4.57", //kDTS_JpegLossless (table A.4-3) "1.2.840.10008.1.2.4.70", //kDTS_JpegLosslessFirstOrder (table A.4-3) "1.2.840.10008.1.2.4.80", //kDTS_JpegLSLossless (A.4.3) "1.2.840.10008.1.2.4.81", //kDTS_JpegLSNearLossless (A.4.3) "1.2.840.10008.1.2.4.90", //kDTS_Jpeg2000LosslessReversible (A.4.4) "1.2.840.10008.1.2.4.91", //kDTS_Jpeg2000LosslessOrLossy (A.4.4) "1.2.840.10008.1.2.4.92", //kDTS_Jpeg2000Part2LosslessReversible (A.4.4) "1.2.840.10008.1.2.4.93", //kDTS_Jpeg2000Part2LosslessOrLossy (A.4.4) "1.2.840.10008.1.2.4.100", //kDTS_MPEG2_ML (table A.4.5) "1.2.840.10008.1.2.4.101", //kDTS_MPEG2_HL (table A.4.5) "1.2.840.10008.1.2.4.102", //kDTS_MPEG4 (table A.4.6) "1.2.840.10008.1.2.4.103", //kDTS_MPEG4Temporal (table A.4.6) "1.2.840.10008.1.2.1.99" //kDTS_Deflate (A.5) //A.6 - DICOM JPIP REFERENCED - 1.2.840.10008.1.2.4.94 //A.7 - DICOM JPIP REFERENCED DEFLATE - 1.2.840.10008.1.2.4.95 }; const EDicomTransferSyntax find_transfer_syntax_for_string(const char *pSyntaxString) { for (int tsIndex = 0; tsIndex < kDTS_Count; ++tsIndex) { if (!strcmp(skTransferSyntaxStrings[tsIndex], pSyntaxString)) { return EDicomTransferSyntax(tsIndex); } } return kDTS_Unknown; } const bool is_jpeg_transfer_syntax(const EDicomTransferSyntax transferSyntax) { return (transferSyntax >= kDTS_JpegBaseline && transferSyntax <= kDTS_Jpeg2000Part2LosslessOrLossy); } //see part 3, C.7.6.3.1.2 enum EDicomPhotometricInterpretation { kDPI_Unknown = 0, kDPI_Monochrome1, kDPI_Monochrome2, kDPI_PaletteColor, kDPI_RGB, kDPI_HSV, //"retired" kDPI_ARGB, //"retired" kDPI_CMYK, //"retired" kDPI_YBR_Full, kDPI_YBR_Full_422, kDPI_YBR_Partial_422, kDPI_YBR_Partial_420, kDPI_YBR_ICT, kDPI_YBR_RCT, kDPI_Count }; static const char *skPhotometricInterpretationStrings[kDPI_Count] = { "Unknown", //kDPI_Unknown "MONOCHROME1", //kDPI_Monochrome1 "MONOCHROME2", //kDPI_Monochrome2 "PALETTE COLOR", //kDPI_PaletteColor "RGB", //kDPI_RGB "HSV", //kDPI_HSV "ARGB", //kDPI_ARGB "CMYK", //kDPI_CMYK "YBR_FULL", //kDPI_YBR_Full "YBR_FULL_422", //kDPI_YBR_Full_422 "YBR_PARTIAL_422", //kDPI_YBR_Partial_422 "YBR_PARTIAL_420", //kDPI_YBR_Partial_420 "YBR_ICT", //kDPI_YBR_ICT "YBR_RCT" //kDPI_YBR_RCT }; const EDicomPhotometricInterpretation photometric_interpretation_for_string(const char *pPhotoString) { for (int piIndex = 0; piIndex < kDPI_Count; ++piIndex) { if (!stricmp(skPhotometricInterpretationStrings[piIndex], pPhotoString)) { return EDicomPhotometricInterpretation(piIndex); } } return kDPI_Unknown; } //see table 6.2-1 for a full list of value representations static const unsigned short skVR_ApplicationEntity = 'EA'; //AE - application entity static const unsigned short skVR_AgeString = 'SA'; //AS - age string static const unsigned short skVR_AttributeTag = 'TA'; //AT - attribute tag static const unsigned short skVR_CodeString = 'SC'; //CS - code string static const unsigned short skVR_Date = 'AD'; //DA - date static const unsigned short skVR_DecimalString = 'SD'; //DS - decimal string static const unsigned short skVR_DateTime = 'TD'; //DT - date time static const unsigned short skVR_Float = 'LF'; //FL - float 32 static const unsigned short skVR_Double = 'DF'; //FD - float 64 static const unsigned short skVR_IntString = 'SI'; //IS - integer string (string character representation of int) static const unsigned short skVR_LongString = 'OL'; //LO - long string (string character representation of long) static const unsigned short skVR_LongText = 'TL'; //LT - long text static const unsigned short skVR_OtherByte = 'BO'; //OB - other byte string static const unsigned short skVR_OtherFloat = 'FO'; //OF - other float string static const unsigned short skVR_OtherWord = 'WO'; //OW - other word (16-bit) string static const unsigned short skVR_PersonName = 'NP'; //PN - person name static const unsigned short skVR_ShortString = 'HS'; //SH - short string (string character representation of short) static const unsigned short skVR_SignedLong = 'LS'; //SL - signed long (32-bit) static const unsigned short skVR_Sequence = 'QS'; //SQ - sequence, see table 7.5-3 for data layout static const unsigned short skVR_SignedShort = 'SS'; //SS - signed short (16-bit) static const unsigned short skVR_ShortText = 'TS'; //ST - short text static const unsigned short skVR_Time = 'MT'; //TM - time static const unsigned short skVR_UID = 'IU'; //UI - unique identifier static const unsigned short skVR_UnsignedLong = 'LU'; //UL - unsigned long (32-bit) static const unsigned short skVR_Unknown = 'NU'; //UN - unknown static const unsigned short skVR_UnsignedShort = 'SU'; //US - unsigned short (16-bit) static const unsigned short skVR_UnlimitedText = 'TU'; //UT - unlimited text static const unsigned short skVR_Invalid = 0; //see section 7.1.2 (table 7.1-1) static const unsigned short skExplicitVRsFor32BitLength[] = { skVR_OtherByte, skVR_OtherWord, skVR_OtherFloat, skVR_Sequence, skVR_UnlimitedText, skVR_Unknown }; static const int skExplicitVRsFor32BitLengthCount = sizeof(skExplicitVRsFor32BitLength) / sizeof(unsigned short); static const int skHeaderOffset = 128; static const int skDataElemsOffset = skHeaderOffset + 4; static const unsigned int skTemporaryStringElementValueSize = 512; static const unsigned short skGroup2 = 0x0002; static const unsigned short skG2TransferSyntaxUID = 0x0010; static const unsigned short skGroupDelimter = 0xFFFE; #define MAKE_ELEMENT_ID(groupNum, elemNum) (groupNum | (elemNum << 16)) static const unsigned int skFrameTimeElement = MAKE_ELEMENT_ID(0x0018, 0x1063); static const unsigned int skWaveformDataElement = MAKE_ELEMENT_ID(0x5400, 0x1010); static const unsigned int skRedPaletteLookupTableDataElement = MAKE_ELEMENT_ID(0x0028, 0x1201); static const unsigned int skGreenPaletteLookupTableDataElement = MAKE_ELEMENT_ID(0x0028, 0x1202); static const unsigned int skBluePaletteLookupTableDataElement = MAKE_ELEMENT_ID(0x0028, 0x1203); static const unsigned int skAlphaPaletteLookupTableDataElement = MAKE_ELEMENT_ID(0x0028, 0x1204); static const unsigned int skRedPaletteLookupTableDescElement = MAKE_ELEMENT_ID(0x0028, 0x1101); static const unsigned int skGreenPaletteLookupTableDescElement = MAKE_ELEMENT_ID(0x0028, 0x1102); static const unsigned int skBluePaletteLookupTableDescElement = MAKE_ELEMENT_ID(0x0028, 0x1103); static const unsigned int skAlphaPaletteLookupTableDescElement = MAKE_ELEMENT_ID(0x0028, 0x1104); static const unsigned int skRedSegPaletteLookupTableDataElement = MAKE_ELEMENT_ID(0x0028, 0x1221); static const unsigned int skGreenSegPaletteLookupTableDataElement = MAKE_ELEMENT_ID(0x0028, 0x1222); static const unsigned int skBlueSegPaletteLookupTableDataElement = MAKE_ELEMENT_ID(0x0028, 0x1223); static const unsigned int skLookupTableDataElement = MAKE_ELEMENT_ID(0x0028, 0x3006); static const unsigned int skLookupTableDescElement = MAKE_ELEMENT_ID(0x0028, 0x3002); static const unsigned int skBlendingLookupTableDataElement = MAKE_ELEMENT_ID(0x0028, 0x1408); static const unsigned int skOverlayDataBaseElement = MAKE_ELEMENT_ID(0x6000, 0x3000); //60xx represents xx overlays //see 8.1.1 //however, pixel macro attributes are only vaguely referenced in part 5, see part 3 table C.7-11b for a complete list. //see also part 3 table C.7-14 for multi-frame attributes. static const unsigned int skPixelDataSamplesPerPixelElement = MAKE_ELEMENT_ID(0x0028, 0x0002); static const unsigned int skPixelDataPhotometricInterpretationElement = MAKE_ELEMENT_ID(0x0028, 0x0004); static const unsigned int skPixelDataRowsElement = MAKE_ELEMENT_ID(0x0028, 0x0010); static const unsigned int skPixelDataColumnsElement = MAKE_ELEMENT_ID(0x0028, 0x0011); static const unsigned int skPixelDataBitsAllocatedElement = MAKE_ELEMENT_ID(0x0028, 0x0100); static const unsigned int skPixelDataBitsStoredElement = MAKE_ELEMENT_ID(0x0028, 0x0101); static const unsigned int skPixelDataHighBitElement = MAKE_ELEMENT_ID(0x0028, 0x0102); static const unsigned int skPixelDataRepresentationElement = MAKE_ELEMENT_ID(0x0028, 0x0103); static const unsigned int skPixelDataElement = MAKE_ELEMENT_ID(0x7FE0, 0x0010); static const unsigned int skPixelDataPlanarConfigurationElement = MAKE_ELEMENT_ID(0x0028, 0x0006); static const unsigned int skPixelDataFrameCountElement = MAKE_ELEMENT_ID(0x0028, 0x0008); static const unsigned int skPixelDataFrameIncrementPointerElement = MAKE_ELEMENT_ID(0x0028, 0x0009); static const unsigned int skPixelDataAspectRatioElement = MAKE_ELEMENT_ID(0x0028, 0x0034); static const unsigned int skPixelDataSmallestValueElement = MAKE_ELEMENT_ID(0x0028, 0x0106); static const unsigned int skPixelDataLargestValueElement = MAKE_ELEMENT_ID(0x0028, 0x0107); static const unsigned int skPixelDataICCProfileElement = MAKE_ELEMENT_ID(0x0028, 0x2000); static const unsigned int skPixelDataRescaleIntercept = MAKE_ELEMENT_ID(0x0028, 0x1052); static const unsigned int skPixelDataRescaleSlope = MAKE_ELEMENT_ID(0x0028, 0x1053); //see part 3 table C.7.6.24-1 (floating point image pixel module attributes) static const unsigned int skFloatPixelDataElement = MAKE_ELEMENT_ID(0x7FE0, 0x0008); static const unsigned int skDoublePixelDataElement = MAKE_ELEMENT_ID(0x7FE0, 0x0009); static const unsigned int skSequenceTag = MAKE_ELEMENT_ID(0xFFFE, 0xE000); static const unsigned int skSequenceEndTag = MAKE_ELEMENT_ID(0xFFFE, 0xE0DD); //possible todo - support volume data, see part 3 table C.8.12.4-1 //not thread-safe const char *temporary_string_for_element_value(RichBitStreamEx &bs, const unsigned int valueLength) { static char sReadBuffer[skTemporaryStringElementValueSize]; const int readLength = std::min<unsigned int>(valueLength, skTemporaryStringElementValueSize-1); bs.ReadBytes(sReadBuffer, readLength); sReadBuffer[readLength] = 0; //eliminate trailing whitespace for (int trailingWhiteSpace = readLength - 1; trailingWhiteSpace > 0; --trailingWhiteSpace) { if (sReadBuffer[trailingWhiteSpace] == ' ' || sReadBuffer[trailingWhiteSpace] == '\r' || sReadBuffer[trailingWhiteSpace] == '\n' || sReadBuffer[trailingWhiteSpace] == 9 /*tab*/) { sReadBuffer[trailingWhiteSpace] = 0; } else if (sReadBuffer[trailingWhiteSpace] != 0) { break; } } return sReadBuffer; } template<class DestType> static const DestType read_single_element_value_for_vr(RichBitStreamEx &bs, const unsigned short vr, const unsigned int valueLength) { switch (vr) { case skVR_Double: return (DestType)bs.ReadDouble(); case skVR_Float: case skVR_OtherFloat: return (DestType)bs.ReadFloat(); case skVR_SignedLong: return (DestType)bs.ReadInt(); case skVR_UnsignedLong: return (DestType)bs.ReadUInt(); case skVR_SignedShort: return (DestType)bs.ReadShort(); case skVR_UnsignedShort: case skVR_OtherWord: case skVR_AttributeTag: return (DestType)bs.ReadUShort(); case skVR_OtherByte: return (DestType)bs.ReadByte(); case skVR_IntString: case skVR_LongString: case skVR_ShortString: { const char *pTempBuffer = temporary_string_for_element_value(bs, valueLength); return (DestType)atoi(pTempBuffer); } case skVR_DecimalString: { const char *pTempBuffer = temporary_string_for_element_value(bs, valueLength); return (DestType)atof(pTempBuffer); } case skVR_Invalid: { //we don't have a VR for this element, try to intuit it based on size switch (valueLength) { case 1: return (DestType)bs.ReadByte(); case 2: return (DestType)bs.ReadUShort(); case 4: return (DestType)bs.ReadUInt(); default: return 0; } } default: return 0; } } //part 3, C.7.6.3.1.5 struct SDicomPalette { SDicomPalette() : mEntryCount(0) , mFirstMappedEntry(0) , mBitsPerEntry(0) , mpData(NULL) , mDataLength(0) , mDataVR(skVR_Invalid) { } union { unsigned int mDescriptor[3]; struct { unsigned int mEntryCount; unsigned int mFirstMappedEntry; unsigned int mBitsPerEntry; }; }; unsigned char *mpData; int mDataLength; unsigned short mDataVR; }; struct SDicomImage { SDicomImage() : mpData(NULL) , mDataLength(0) , mDataVR(skVR_Invalid) , mSamplesPerPixel(0) , mBitsPerChannel(0) , mBitsStoredPerChannel(0) , mHighBit(0) , mPhotometricInterpretation(kDPI_Unknown) , mPlanarConfiguration(0) , mRowCount(0) , mColumnCount(0) , mSliceCount(0) , mRescaleIntercept(0.0f) , mRescaleSlope(1.0f) { } unsigned char *mpData; int mDataLength; //we should always end up with an explicit VR for image data, even with an implicit VR syntax, //because implicit VR requirements dictate OW. unsigned short mDataVR; int mSamplesPerPixel; int mBitsPerChannel; int mBitsStoredPerChannel; int mHighBit; EDicomPhotometricInterpretation mPhotometricInterpretation; int mPlanarConfiguration; int mRowCount; int mColumnCount; int mSliceCount; //scale and bias on data (intercept is bias, slope is scale) float mRescaleIntercept; float mRescaleSlope; SDicomPalette mPaletteRgba[4]; }; class CDicomDataElement { public: //transfer syntax may be provided as kDTS_Unknown before a skG2TransferSyntaxUID element is found. explicit CDicomDataElement(RichBitStreamEx &bs, const EDicomTransferSyntax transferSyntax) { bs.SetFlags(bs.GetFlags() & ~BITSTREAMFL_BIGENDIAN); mGroupNum = bs.ReadUShort(); //set the stream endian for the rest of the element header and data if (SetStreamEndian(bs, transferSyntax) && mGroupNum != skGroup2) { //group numbers other than 2 are endian-swapped LITTLE_BIG_SWAP(mGroupNum); } mElemNum = bs.ReadUShort(); bool isStandardElem = true; const bool isImplicitVR = (transferSyntax == kDTS_ImplicitVRLittleEndian && mGroupNum != skGroup2); if (isImplicitVR) { //see Annex A (A.1) which specifies implicit VR requirements for these groups/elements. //many of these same requirements are dictated for default explicit VR, but we don't really //give a shit if the binary wants to specify something non-standards-conformant. if ((mElementId & 0xFFFFFF00) == skOverlayDataBaseElement) { //special-case for overlay data mVR = skVR_OtherWord; } else { switch (mElementId) { case skPixelDataElement: case skWaveformDataElement: case skRedPaletteLookupTableDataElement: case skGreenPaletteLookupTableDataElement: case skBluePaletteLookupTableDataElement: case skAlphaPaletteLookupTableDataElement: case skRedSegPaletteLookupTableDataElement: case skGreenSegPaletteLookupTableDataElement: case skBlueSegPaletteLookupTableDataElement: case skLookupTableDataElement: case skBlendingLookupTableDataElement: mVR = skVR_OtherWord; break; case skRedPaletteLookupTableDescElement: case skGreenPaletteLookupTableDescElement: case skBluePaletteLookupTableDescElement: case skAlphaPaletteLookupTableDescElement: case skLookupTableDescElement: mVR = skVR_UnsignedShort; break; case skPixelDataRescaleIntercept: case skPixelDataRescaleSlope: mVR = skVR_DecimalString; break; default: mVR = skVR_Invalid; break; } } } else if (mGroupNum == skGroupDelimter) { mVR = skVR_Invalid; } else { //read VR byte by byte as it's actually a string representation mVR = bs.ReadByte() | ((unsigned short)bs.ReadByte() << 8); isStandardElem = !ExplicitVRFor32BitLength(); } if (isStandardElem) { //see 7.1.3, length is always 32 bits in implicit VR mode and for delimiters mValueLength = (isImplicitVR || mGroupNum == skGroupDelimter) ? bs.ReadUInt() : bs.ReadUShort(); } else { bs.ReadUShort(); //reserved bytes, discard mValueLength = bs.ReadUInt(); } mValueOfs = bs.GetOffset(); if (mGroupNum == skGroupDelimter || mVR == skVR_Sequence) { //we don't care about respecting sequences encountered mid-stream, so just skip over the //element header. mValueLength = 0; } else if (mValueLength == 0xFFFFFFFF) { //section 7.1.1 notes that SQ and UN can be total shitfaces about providing a value length. //7.5 has further details on item encoding and delimitation. //it seems there are also some proprietary data element types which can do this. OW and OB //may also be "undefined" depending on transfer syntax. if this happens, we try to iterate //over the sequence tags to find the total length. if we don't start on a sequence tag, we //just say fuck it and act like the start of the data is the start of the next data element, //which may or may not explode gloriously under unexpected circumstances. unsigned int delimiter = bs.ReadUInt(); if (delimiter != skSequenceTag) { //alright then, fuck it. mValueLength = 0; } else { while (delimiter == skSequenceTag || delimiter == skSequenceEndTag) { const int seqSegSize = bs.ReadInt(); if (seqSegSize < 0) { break; } bs.SetOffset(bs.GetOffset() + seqSegSize); if (delimiter == skSequenceEndTag) { //all done break; } delimiter = bs.ReadUInt(); } mValueLength = bs.GetOffset() - mValueOfs; } //default us back to the start of the value data bs.SetOffset(mValueOfs); } } explicit CDicomDataElement(const unsigned int elementId, const unsigned short vr) : mElementId(elementId) , mVR(vr) , mValueLength(0) , mValueOfs(0) { } explicit CDicomDataElement(const unsigned short groupNum, const unsigned short elemNum, const unsigned short vr) : mGroupNum(groupNum) , mElemNum(elemNum) , mVR(vr) , mValueLength(0) , mValueOfs(0) { } void PrepDataStreamForValueRead(RichBitStreamEx &bs, const EDicomTransferSyntax transferSyntax) const { //set the correct endian mode given the group number and uid-based preference SetStreamEndian(bs, transferSyntax); bs.SetOffset(mValueOfs); } void SetStreamToNextDataElement(RichBitStreamEx &bs) const { bs.SetOffset(mValueOfs + mValueLength); } template<class DestType> const DestType ReadSingleElementValue(RichBitStreamEx &bs, const EDicomTransferSyntax transferSyntax) const { PrepDataStreamForValueRead(bs, transferSyntax); return read_single_element_value_for_vr<DestType>(bs, mVR, mValueLength); } template<class DestType> void ReadArrayElementValues(DestType *pOutValues, const int valueCount, RichBitStreamEx &bs, const EDicomTransferSyntax transferSyntax) const { PrepDataStreamForValueRead(bs, transferSyntax); for (int valueIndex = 0; valueIndex < valueCount; ++valueIndex) { pOutValues[valueIndex] = read_single_element_value_for_vr<DestType>(bs, mVR, mValueLength); } } //not thread-safe const char *TemporaryStringForElementValue(RichBitStreamEx &bs) const { return temporary_string_for_element_value(bs, mValueLength); } const unsigned short GetGroupNum() const { return mGroupNum; } const unsigned short GetElemNum() const { return mElemNum; } const unsigned int GetElementId() const { return mElementId; } const unsigned short GetVR() const { return mVR; } const unsigned int GetValueLength() const { return mValueLength; } const unsigned int GetValueOfs() const { return mValueOfs; } void WriteToStream(RichBitStreamEx &bs, const EDicomTransferSyntax transferSyntax, const void *pData, const int dataSize) const { //possible todo - do the right thing here based on transfer syntax bs.WriteUShort(mGroupNum); bs.WriteUShort(mElemNum); const unsigned char *pVR = (const unsigned char *)&mVR; bs.WriteByte(pVR[0]); bs.WriteByte(pVR[1]); const int alignedSize = ((dataSize + 1) & ~1); if (!ExplicitVRFor32BitLength()) { bs.WriteUShort((unsigned short)alignedSize); } else { bs.WriteUShort(0); //reserved bytes bs.WriteUInt(alignedSize); } if (pData) { //write the data bs.WriteBytes(pData, dataSize); if (alignedSize > dataSize) { //write even byte pad bs.WriteByte(0); } } } template<class DataType> void WriteTypeToStream(RichBitStreamEx &bs, const EDicomTransferSyntax transferSyntax, const DataType &data) const { WriteToStream(bs, transferSyntax, &data, sizeof(DataType)); } void WriteStringToStream(RichBitStreamEx &bs, const EDicomTransferSyntax transferSyntax, const char *pDataString) const { //terminator not included const int stringLength = (int)strlen(pDataString); WriteToStream(bs, transferSyntax, pDataString, stringLength); } private: bool SetStreamEndian(RichBitStreamEx &bs, const EDicomTransferSyntax transferSyntax) const { bool isBigEndian = (transferSyntax == kDTS_ExplicitVRBigEndian); const int flags = bs.GetFlags(); const int streamEndianFlags = (mGroupNum == skGroup2 || !isBigEndian) ? (flags & ~BITSTREAMFL_BIGENDIAN) : (flags | BITSTREAMFL_BIGENDIAN); bs.SetFlags(streamEndianFlags); return isBigEndian; } bool ExplicitVRFor32BitLength() const { for (int vrCheckIndex = 0; vrCheckIndex < skExplicitVRsFor32BitLengthCount; ++vrCheckIndex) { if (skExplicitVRsFor32BitLength[vrCheckIndex] == mVR) { //value representation is in the list of types with reserved bytes and 32-bit length return true; } } return false; } //see figure 7.1-1 union { unsigned int mElementId; struct { unsigned short mGroupNum; unsigned short mElemNum; }; }; unsigned short mVR; unsigned int mValueLength; unsigned int mValueOfs; //offset from the start of stream data, in bytes }; static const EDicomTransferSyntax find_transfer_syntax(RichBitStreamEx &bs) { while (bs.GetOffset() < bs.GetSize()) { //iterate the group2 elements until we find the transfer syntax const CDicomDataElement elem(bs, kDTS_Unknown); if (elem.GetGroupNum() != skGroup2) { //always expect group2 elements up to the uid return kDTS_Unknown; } const int endOfElemData = int(elem.GetValueOfs() + elem.GetValueLength()); if (endOfElemData <= 0 || endOfElemData >= bs.GetSize()) { //bad data size return kDTS_Unknown; } if (elem.GetElemNum() == skG2TransferSyntaxUID) { //alright, found it. let's see what we're dealing with. bs.SetOffset(elem.GetValueOfs()); const char *pTransferSyntaxString = elem.TemporaryStringForElementValue(bs); //leave the stream at the start of the next element elem.SetStreamToNextDataElement(bs); return find_transfer_syntax_for_string(pTransferSyntaxString); } else { //skip it elem.SetStreamToNextDataElement(bs); } } return kDTS_Unknown; } bool Image_MRI_CheckDicom(BYTE *fileBuffer, int bufferLen, noeRAPI_t *rapi) { if (bufferLen <= skDataElemsOffset) { return false; } if (memcmp(fileBuffer + skHeaderOffset, "DICM", 4) != 0) { return false; } RichBitStreamEx bs(fileBuffer, bufferLen); bs.SetOffset(skDataElemsOffset); const EDicomTransferSyntax transferSyntax = find_transfer_syntax(bs); return transferSyntax != kDTS_Unknown; } static void converted_float_to_rgb(unsigned char *pRgbOut, const int channelsOut, const int width, const int height, const float *pConvertedFloat, const int channelsPerPixel, const unsigned short sourceVR, const float minSampleVal, const float maxSampleVal) { const float sampleRange = (maxSampleVal - minSampleVal); for (int pixelIndex = 0; pixelIndex < width * height; ++pixelIndex) { unsigned char *pRgbDst = pRgbOut + pixelIndex * channelsOut; const float *pFloatSrc = pConvertedFloat + pixelIndex * channelsPerPixel; for (int channelIndex = 0; channelIndex < channelsPerPixel; ++channelIndex) { if (sourceVR == skVR_OtherByte) { //if the source was u8 data, convert 1:1 pRgbDst[channelIndex] = (unsigned char)pFloatSrc[channelIndex]; } else { //otherwise, quantize into the utilized range const float v = (pFloatSrc[channelIndex] - minSampleVal) / sampleRange; pRgbDst[channelIndex] = (unsigned char)(v * 255.0f); } } //copy up from 0 if the float source doesn't have enough channels for (int remainingChannelIndex = channelsPerPixel; remainingChannelIndex < channelsOut; ++remainingChannelIndex) { pRgbDst[remainingChannelIndex] = (remainingChannelIndex == 3) ? 255 : pRgbDst[0]; } } } static float normalize_value_for_image_data(const float value, const SDicomImage &dicomImage) { float valueOut = value; if (dicomImage.mDataVR == skVR_OtherWord && dicomImage.mBitsPerChannel > 0) { valueOut /= (float)((1 << dicomImage.mBitsPerChannel) - 1); } else { float scale; switch (dicomImage.mDataVR) { case skVR_SignedLong: scale = 1.0f / 2147483647.0f; break; case skVR_UnsignedLong: scale = 1.0f / 4294967295.0f; break; case skVR_SignedShort: scale = 1.0f / 32767.0f; break; case skVR_UnsignedShort: case skVR_OtherWord: scale = 1.0f / 65535.0f; break; case skVR_OtherByte: scale = 1.0f / 255.0f; break; default: scale = 1.0f; break; } valueOut *= scale; } return valueOut; } static float apply_value_slope(const float value, const SDicomImage &dicomImage) { return (value * dicomImage.mRescaleSlope) + dicomImage.mRescaleIntercept; } static float apply_value_sbp(const float value) { //apply bias first, as we typically want to use it to range-correct return powf((g_dicomOpts->mSBPBias + value) * g_dicomOpts->mSBPScale, g_dicomOpts->mSBPExponent); } static void convert_dicom_image(CArrayList<noesisTex_t *> &texturesOut, const SDicomImage &dicomImage, const EDicomTransferSyntax transferSyntax, noeRAPI_t *rapi) { if (!dicomImage.mpData || dicomImage.mDataLength <= 0 || dicomImage.mDataVR == skVR_Invalid || dicomImage.mBitsPerChannel <= 0) { rapi->LogOutput("WARNING: Invalid image data/size, discarding.\n"); return; } if (dicomImage.mRowCount <= 0 || dicomImage.mColumnCount <= 0) { rapi->LogOutput("WARNING: Invalid image dimensions, discarding.\n"); return; } //don't do anything intelligent about orientation at the moment const int width = dicomImage.mColumnCount; const int height = dicomImage.mRowCount; const int sliceCount = std::max<int>(dicomImage.mSliceCount, 1); const bool applySlope = (g_dicomOpts && g_dicomOpts->mApplySlope); const bool applySBP = (g_dicomOpts && g_dicomOpts->mApplySBP); const bool dataNormalize = (g_dicomOpts && g_dicomOpts->mDataNormalize); RichBitStreamEx bs(dicomImage.mpData, dicomImage.mDataLength); //possible todo - do we need to switch the endianness of the image data stream with a big-endian transfer syntax? //writing this to parse each pixel element through the bitstream in case this does end up being necessary. for (int sliceIndex = 0; sliceIndex < sliceCount; ++sliceIndex) { unsigned char *pConvertedRgb = NULL; float *pConvertedFloat = NULL; ENoeHdrTexFormat convertedFloatFormat; bool convertedHasAlpha = false; bool isColor = false; bool expectJpeg = false; switch (dicomImage.mPhotometricInterpretation) { //we rely on libjpeg to do the color conversion for these interpretations, and just accept them here as rgb. //if we came across one in non-jpeg form, this would certainly not work! case kDPI_YBR_Full: case kDPI_YBR_Full_422: case kDPI_YBR_Partial_422: case kDPI_YBR_Partial_420: case kDPI_YBR_ICT: case kDPI_YBR_RCT: expectJpeg = true; //fall through intentionally case kDPI_RGB: isColor = true; //fall through intentionally case kDPI_Monochrome1: case kDPI_Monochrome2: { if (expectJpeg && !is_jpeg_transfer_syntax(transferSyntax)) { rapi->LogOutput("WARNING: Photometric interpretation not supported for transfer syntax.\n"); break; } const int channelsPerPixel = (isColor) ? 3 : 1; convertedFloatFormat = (isColor) ? kNHDRTF_RGB_F96 : kNHDRTF_Lum_F32; pConvertedFloat = (float *)rapi->Noesis_UnpooledAlloc(width * height * sizeof(float) * channelsPerPixel); pConvertedRgb = (unsigned char *)rapi->Noesis_UnpooledAlloc(width * height * 3); float minSampleVal = FLT_MAX; float maxSampleVal = -FLT_MAX; //possible todo - could optimize this, written this way for easy tinkering const bool planarMode = (dicomImage.mPlanarConfiguration != 0); const int channelsOuterLoop = (planarMode) ? channelsPerPixel : 1; const int channelsInnerLoop = (!planarMode) ? channelsPerPixel : 1; for (int channelIndexOuter = 0; channelIndexOuter < channelsOuterLoop; ++channelIndexOuter) { for (int y = 0; y < height; ++y) { for (int x = 0; x < width; ++x) { const int pixelIndex = y * width + x; for (int channelIndexInner = 0; channelIndexInner < channelsInnerLoop; ++channelIndexInner) { float &destFloat = pConvertedFloat[pixelIndex * channelsPerPixel + channelIndexInner + channelIndexOuter]; if (dicomImage.mDataVR == skVR_OtherWord && dicomImage.mBitsPerChannel > 0) { //do an explicit bit read destFloat = (float)bs.ReadBits(dicomImage.mBitsPerChannel); } else { destFloat = read_single_element_value_for_vr<float>(bs, dicomImage.mDataVR, 0); } if (dataNormalize) { destFloat = normalize_value_for_image_data(destFloat, dicomImage); } if (applySlope) { destFloat = apply_value_slope(destFloat, dicomImage); } if (applySBP) { destFloat = apply_value_sbp(destFloat); } minSampleVal = (destFloat < minSampleVal) ? destFloat : minSampleVal; maxSampleVal = (destFloat > maxSampleVal) ? destFloat : maxSampleVal; } } } } //sanity-check range if (minSampleVal >= maxSampleVal) { minSampleVal = 0.0f; maxSampleVal = 1.0f; } converted_float_to_rgb(pConvertedRgb, 3, width, height, pConvertedFloat, channelsPerPixel, dicomImage.mDataVR, minSampleVal, maxSampleVal); } break; case kDPI_PaletteColor: { const SDicomPalette *pPalettes = dicomImage.mPaletteRgba; if (pPalettes[0].mpData && pPalettes[0].mBitsPerEntry > 0 && pPalettes[0].mEntryCount > 0 && pPalettes[1].mpData && pPalettes[1].mBitsPerEntry > 0 && pPalettes[1].mEntryCount > 0 && pPalettes[2].mpData && pPalettes[2].mBitsPerEntry > 0 && pPalettes[2].mEntryCount > 0) { const bool alphaIsValid = (pPalettes[3].mpData && pPalettes[3].mBitsPerEntry > 0 && pPalettes[3].mEntryCount > 0); float *pPaletteColors[4] = { (float *)rapi->Noesis_UnpooledAlloc(sizeof(float) * pPalettes[0].mEntryCount), (float *)rapi->Noesis_UnpooledAlloc(sizeof(float) * pPalettes[1].mEntryCount), (float *)rapi->Noesis_UnpooledAlloc(sizeof(float) * pPalettes[2].mEntryCount), alphaIsValid ? (float *)rapi->Noesis_UnpooledAlloc(sizeof(float) * pPalettes[3].mEntryCount) : NULL }; //convert palette colors for (int paletteIndex = 0; paletteIndex < 4; ++paletteIndex) { float *pPaletteDst = pPaletteColors[paletteIndex]; if (pPaletteDst) { const SDicomPalette *pPalette = pPalettes + paletteIndex; //palette colors are always normalized to the available range const float colorScale = (pPalette->mBitsPerEntry >= 16) ? 1.0f / 65535.0f : 1.0f / 255.0f; RichBitStreamEx palBs(pPalette->mpData, pPalette->mDataLength); for (unsigned int entryIndex = 0; entryIndex < pPalette->mEntryCount; ++entryIndex) { float &destFloat = pPaletteDst[entryIndex]; destFloat = read_single_element_value_for_vr<float>(palBs, pPalette->mDataVR, 0) * colorScale; //apply value processing directly to palette values if (applySlope) { destFloat = apply_value_slope(destFloat, dicomImage); } if (applySBP) { destFloat = apply_value_sbp(destFloat); } } } } convertedFloatFormat = kNHDRTF_RGBA_F128; convertedHasAlpha = true; pConvertedFloat = (float *)rapi->Noesis_UnpooledAlloc(width * height * sizeof(float) * 4); pConvertedRgb = (unsigned char *)rapi->Noesis_UnpooledAlloc(width * height * 4); //seems like we ignore the actual VR for the pixel data for this mode, but doing this //check just in case. const bool ushortPixels = (dicomImage.mDataLength == (width * height * 2)); for (int y = 0; y < height; ++y) { for (int x = 0; x < width; ++x) { float *pDestPixel = pConvertedFloat + (y * width + x) * 4; const int palEntryIndex = (dicomImage.mBitsPerChannel > 0) ? bs.ReadBits(dicomImage.mBitsPerChannel) : (ushortPixels) ? bs.ReadUShort() : bs.ReadByte(); for (int paletteIndex = 0; paletteIndex < 4; ++paletteIndex) { const float *pPaletteFloat = pPaletteColors[paletteIndex]; if (pPaletteFloat) { const SDicomPalette *pPalette = pPalettes + paletteIndex; int clampedPalEntryIndex = std::max<int>(palEntryIndex, pPalette->mFirstMappedEntry); clampedPalEntryIndex = std::min<int>(clampedPalEntryIndex, pPalette->mEntryCount - 1); pDestPixel[paletteIndex] = pPaletteFloat[clampedPalEntryIndex]; } else { pDestPixel[paletteIndex] = 1.0f; } } } } rapi->Noesis_UnpooledFree(pPaletteColors[0]); rapi->Noesis_UnpooledFree(pPaletteColors[1]); rapi->Noesis_UnpooledFree(pPaletteColors[2]); if (pPaletteColors[3]) { rapi->Noesis_UnpooledFree(pPaletteColors[3]); } //colors are already in 0..1 const float minSampleVal = 0.0f; const float maxSampleVal = 1.0f; converted_float_to_rgb(pConvertedRgb, 4, width, height, pConvertedFloat, 4, pPalettes->mDataVR, minSampleVal, maxSampleVal); } else { rapi->LogOutput("WARNING: Photometric interpretation is %s, but we're missing palette data.\n", skPhotometricInterpretationStrings[dicomImage.mPhotometricInterpretation]); } } break; default: rapi->LogOutput("WARNING: Photometric interpretation not currently supported: %s\n", skPhotometricInterpretationStrings[dicomImage.mPhotometricInterpretation]); break; } if (!pConvertedRgb) { break; } char sliceName[512]; sprintf_s(sliceName, "image%04i", texturesOut.Num()); const int colorChannelCount = (convertedHasAlpha) ? 4 : 3; noesisTex_t *pTex = rapi->Noesis_TextureAllocEx(sliceName, width, height, pConvertedRgb, width * height * colorChannelCount, (convertedHasAlpha) ? NOESISTEX_RGBA32 : NOESISTEX_RGB24, 0, 1); if (pConvertedFloat) { int floatChannelCount; switch (convertedFloatFormat) { case kNHDRTF_RGB_F96: floatChannelCount = 3; break; case kNHDRTF_RGBA_F128: floatChannelCount = 4; break; case kNHDRTF_Lum_F32: default: floatChannelCount = 1; break; } if (g_dicomOpts && g_dicomOpts->mFloatNormalize) { //normalize floating point data if desired float maxValue = 0.0f; for (int valueIndex = 0; valueIndex < width * height * floatChannelCount; ++valueIndex) { const float v = pConvertedFloat[valueIndex]; maxValue = (v > maxValue) ? v : maxValue; } if (maxValue > 0.0f) { const float invMaxValue = 1.0f / maxValue; for (int valueIndex = 0; valueIndex < width * height * floatChannelCount; ++valueIndex) { pConvertedFloat[valueIndex] *= invMaxValue; } } } pTex->pHdr = rapi->Noesis_AllocHDRTexStructure(pConvertedFloat, width * height * sizeof(float) * floatChannelCount, convertedFloatFormat, NULL, NULL); } pTex->shouldFreeData = true; texturesOut.Append(pTex); } } static unsigned short *jpeg_decode_ls(const void *pSource, const int sourceSize, const int dataPrecision, const SDicomImage &dicomImage, noeRAPI_t *rapi, int *pWidthOut, int *pHeightOut, int *pComponentsOut, int *pPrecisionOut) { JlsParameters params; if (JpegLsReadHeader(pSource, sourceSize, ¶ms) != OK) { return NULL; } const int destSize = params.bytesperline * params.height; unsigned short *pDest = (unsigned short *)rapi->Noesis_UnpooledAlloc(destSize); //modify parameters for the decode params.bytesperline = sizeof(unsigned short) * params.width; params.bitspersample = 16; if (JpegLsDecode(pDest, destSize, pSource, sourceSize, ¶ms) != OK) { rapi->Noesis_UnpooledFree(pDest); return NULL; } *pWidthOut = params.width; *pHeightOut = params.height; *pComponentsOut = params.components; *pPrecisionOut = params.bitspersample; return pDest; } static void decode_jpeg_to_stream(RichBitStreamEx &bsJpegConcat, RichBitStreamEx &bsTotalOut, const SDicomImage &dicomImage, const EDicomTransferSyntax transferSyntax, noeRAPI_t *rapi) { const int jpegSize = bsJpegConcat.GetOffset(); if (jpegSize > 0) { int width, height, components, decodedPrecision; const int dataPrecision = (dicomImage.mBitsStoredPerChannel > 0) ? dicomImage.mBitsStoredPerChannel : 8; unsigned short *pDecoded; switch (transferSyntax) { case kDTS_JpegLSLossless: case kDTS_JpegLSNearLossless: pDecoded = jpeg_decode_ls((const unsigned char *)bsJpegConcat.GetBuffer(), jpegSize, dataPrecision, dicomImage, rapi, &width, &height, &components, &decodedPrecision); break; case kDTS_Jpeg2000LosslessReversible: case kDTS_Jpeg2000LosslessOrLossy: case kDTS_Jpeg2000Part2LosslessReversible: case kDTS_Jpeg2000Part2LosslessOrLossy: pDecoded = Image_JPEG2000_ReadDirectU16((const unsigned char *)bsJpegConcat.GetBuffer(), jpegSize, &width, &height, &components, &decodedPrecision, rapi); break; default: pDecoded = rapi->Image_JPEG_ReadDirect((const unsigned char *)bsJpegConcat.GetBuffer(), jpegSize, dataPrecision, &width, &height, &components, &decodedPrecision, NULL); break; } if (pDecoded) { if (dataPrecision != decodedPrecision) { rapi->LogOutput("WARNING: DICOM stored bits per channel != JPEG data precision: %i vs %i\n", dataPrecision, decodedPrecision); } const int pixelCount = width * height * components; //shift to the dicom-desired range. const int shiftValue = (dicomImage.mBitsPerChannel - decodedPrecision); for (int pixelIndex = 0; pixelIndex < pixelCount; ++pixelIndex) { int pixelValue = pDecoded[pixelIndex]; if (shiftValue > 0) { pixelValue <<= shiftValue; } else { pixelValue >>= (-shiftValue); } bsTotalOut.WriteBits(pixelValue, dicomImage.mBitsPerChannel); } rapi->Noesis_UnpooledFree(pDecoded); } bsJpegConcat.SetOffset(0); } } static unsigned char *decompress_jpeg(int *pSizeOut, const unsigned char *pData, const SDicomImage &dicomImage, const EDicomTransferSyntax transferSyntax, const int maxDataSize, noeRAPI_t *rapi) { //hacked together bullshit, can't find any decent documentation on this. RichBitStreamEx bsTotalOut; RichBitStreamEx bsJpegConcat; RichBitStreamEx bs((void *)pData, maxDataSize); unsigned int delimiter = bs.ReadUInt(); bool firstSeg = true; const int *pFrameOffsets = NULL; int frameCount = 0; int currentFrame = 0; int baseOfs = 0; while (delimiter == skSequenceTag && bs.GetOffset() < bs.GetSize()) { const int seqSegSize = bs.ReadInt(); const int jpegHeaderOfs = bs.GetOffset(); if (seqSegSize > 0) { if (firstSeg) { pFrameOffsets = (const int *)((const char *)bs.GetBuffer() + jpegHeaderOfs); frameCount = seqSegSize / sizeof(int); baseOfs = jpegHeaderOfs + seqSegSize + 8; } else { if (frameCount > 0 && currentFrame < frameCount) { //see if we're at the start of a new frame. if so, decode and reset the offset. int currentOfs = jpegHeaderOfs - baseOfs; if (pFrameOffsets[currentFrame] >= currentOfs) { //flush to the main stream decode_jpeg_to_stream(bsJpegConcat, bsTotalOut, dicomImage, transferSyntax, rapi); ++currentFrame; } } const unsigned char *pJpegData = (const unsigned char *)bs.GetBuffer() + jpegHeaderOfs; bsJpegConcat.WriteBytes(pJpegData, seqSegSize); } } firstSeg = false; bs.SetOffset(jpegHeaderOfs + seqSegSize); delimiter = bs.ReadUInt(); } //flush the rest of the stream decode_jpeg_to_stream(bsJpegConcat, bsTotalOut, dicomImage, transferSyntax, rapi); const int outSize = bsTotalOut.GetOffset(); if (outSize <= 0) { return NULL; } *pSizeOut = outSize; unsigned char *pOutBuffer = (unsigned char *)rapi->Noesis_UnpooledAlloc(outSize); memcpy(pOutBuffer, bsTotalOut.GetBuffer(), outSize); return pOutBuffer; } //see Annex G static unsigned char *decompress_rle(int *pSizeOut, const unsigned char *pData, const int maxDataSize, const int bitsPerChannel, noeRAPI_t *rapi) { RichBitStreamEx bsOut; RichBitStreamEx bs((void *)pData, maxDataSize); unsigned int delimiter = bs.ReadUInt(); while (delimiter == skSequenceTag && bs.GetOffset() < bs.GetSize()) { const unsigned int seqSegSize = bs.ReadUInt(); const int rleHeaderOfs = bs.GetOffset(); const int dataOutOfs = bsOut.GetOffset(); int segCount = bs.ReadInt(); const unsigned int *pSegOffsets = (const unsigned int *)((unsigned char *)bs.GetBuffer() + bs.GetOffset()); for (int segIndex = 0; segIndex < segCount; ++segIndex) { const int segEnd = (segIndex >= (segCount - 1) || pSegOffsets[segIndex + 1] < pSegOffsets[segIndex]) ? rleHeaderOfs + seqSegSize : rleHeaderOfs + pSegOffsets[segIndex + 1]; bs.SetOffset(rleHeaderOfs + pSegOffsets[segIndex]); while (bs.GetOffset() < segEnd) { const char n = bs.ReadChar(); if (n >= 0) { const int repCount = n + 1; for (int valueIndex = 0; valueIndex < repCount && bs.GetOffset() < segEnd; ++valueIndex) { bsOut.WriteByte(bs.ReadByte()); } } else if (n >= -127 && bs.GetOffset() < segEnd) { const unsigned char rep = bs.ReadByte(); const int repCount = -n + 1; for (int valueIndex = 0; valueIndex < repCount; ++valueIndex) { bsOut.WriteByte(rep); } } } } if (segCount > 1) { //now we've gotta interleave the bytes const unsigned int flatDataSize = bsOut.GetOffset() - dataOutOfs; unsigned char *pFlatData = (unsigned char *)bsOut.GetBuffer() + dataOutOfs; unsigned char *pTempBuffer = (unsigned char *)rapi->Noesis_UnpooledAlloc(flatDataSize + segCount); unsigned char *pTempChannels = pTempBuffer + flatDataSize; const int interleavedPixelCount = flatDataSize / segCount; //possible todo - do we need to do bitsPerChannel + 7 here? can we have 9..15-bit rle images? const int bytesPerChannel = bitsPerChannel / 8; //possible todo - optimize this, needs to be refactored after incrementally discovering what a //shitmess the rle is. for (int pixelIndex = 0; pixelIndex < interleavedPixelCount; ++pixelIndex) { //"channel index" isn't necessarily the case, but it probably is. //monochrome 16 bit formats for example will still split each 8 bit //part of each pixel across rle streams. for (int channelIndex = 0; channelIndex < segCount; ++channelIndex) { pTempBuffer[pixelIndex * segCount + channelIndex] = pFlatData[pixelIndex + channelIndex * interleavedPixelCount]; } //apparently we also need to swap bytes back within channels if those //channels are multi-byte. what a beautiful fucking mess! if (bytesPerChannel > 1) { unsigned char *pInterleavedPixel = pTempBuffer + pixelIndex * segCount; for (int channelIndex = 0; channelIndex < segCount; channelIndex += bytesPerChannel) { for (int channelByteIndex = 0; channelByteIndex < bytesPerChannel; ++channelByteIndex) { pTempChannels[channelIndex + channelByteIndex] = pInterleavedPixel[channelIndex + (bytesPerChannel - channelByteIndex - 1)]; } } memcpy(pInterleavedPixel, pTempChannels, segCount); } } //stomp back over the stream data memcpy(pFlatData, pTempBuffer, flatDataSize); rapi->Noesis_UnpooledFree(pTempBuffer); } bs.SetOffset(rleHeaderOfs + seqSegSize); delimiter = bs.ReadUInt(); } const int outSize = bsOut.GetOffset(); if (outSize <= 0) { return NULL; } *pSizeOut = outSize; unsigned char *pOutBuffer = (unsigned char *)rapi->Noesis_UnpooledAlloc(outSize); memcpy(pOutBuffer, bsOut.GetBuffer(), outSize); return pOutBuffer; } static void set_palette_data_from_element(SDicomPalette *pPalette, RichBitStreamEx &bs, const CDicomDataElement &elem) { pPalette->mpData = (unsigned char *)bs.GetBuffer() + elem.GetValueOfs(); pPalette->mDataLength = elem.GetValueLength(); pPalette->mDataVR = elem.GetVR(); } static void set_palette_descriptor_from_element(SDicomPalette *pPalette, RichBitStreamEx &bs, const EDicomTransferSyntax transferSyntax, const CDicomDataElement &elem) { elem.ReadArrayElementValues<unsigned int>(pPalette->mDescriptor, 3, bs, transferSyntax); if (pPalette->mEntryCount == 0) { //special-case - 0 means 2^16 pPalette->mEntryCount = 65536; } } bool Image_MRI_LoadDicom(BYTE *fileBuffer, int bufferLen, CArrayList<noesisTex_t *> &noeTex, noeRAPI_t *rapi) { RichBitStreamEx bs(fileBuffer, bufferLen); bs.SetOffset(skDataElemsOffset); bool transferSyntaxSupported; const EDicomTransferSyntax transferSyntax = find_transfer_syntax(bs); switch (transferSyntax) { case kDTS_ImplicitVRLittleEndian: case kDTS_ExplicitVRLittleEndian: case kDTS_ExplicitVRBigEndian: case kDTS_RLE: case kDTS_JpegBaseline: case kDTS_JpegExtended: case kDTS_JpegSpectral: case kDTS_JpegProgressive: case kDTS_JpegLossless: case kDTS_JpegLosslessFirstOrder: case kDTS_JpegLSLossless: case kDTS_JpegLSNearLossless: case kDTS_Jpeg2000LosslessReversible: case kDTS_Jpeg2000LosslessOrLossy: case kDTS_Jpeg2000Part2LosslessReversible: case kDTS_Jpeg2000Part2LosslessOrLossy: transferSyntaxSupported = true; break; default: transferSyntaxSupported = false; break; } if (!transferSyntaxSupported) { rapi->LogOutput("The following DICOM transfer syntax is not currently supported: %s\n", skTransferSyntaxStrings[transferSyntax]); return NULL; } rapi->LogOutput("Parsing data elements with transfer syntax %s.\n", skTransferSyntaxStrings[transferSyntax]); const int logDataElements = (g_dicomOpts) ? g_dicomOpts->mLogDataElements : 0; if (logDataElements > 0) { //if we want to log all of the data elements, set us back so we iterate all of the g2 elements too. bs.SetOffset(skDataElemsOffset); } int dataElementCount = 0; SDicomImage currentImage; //resume iterating through data elements after the transfer syntax element while (bs.GetOffset() < bs.GetSize()) { const CDicomDataElement elem(bs, transferSyntax); if (logDataElements > 0) { const unsigned short vr = elem.GetVR(); const char *pVR = (vr != 0) ? (const char *)&vr : "NA"; rapi->LogOutput("%04i - Data Element (%04x,%04x) - VR: %c%c Size: %i Offset: %i\n", dataElementCount, elem.GetGroupNum(), elem.GetElemNum(), pVR[0], pVR[1], elem.GetValueLength(), elem.GetValueOfs()); if (logDataElements > 1) { //give data info based on the value representation rapi->LogOutput(" Value: "); switch (elem.GetVR()) { case skVR_ApplicationEntity: case skVR_AgeString: case skVR_CodeString: case skVR_Date: case skVR_DecimalString: case skVR_DateTime: case skVR_IntString: case skVR_LongString: case skVR_LongText: case skVR_PersonName: case skVR_ShortString: case skVR_ShortText: case skVR_Time: case skVR_UID: { const char *pTempBuffer = elem.TemporaryStringForElementValue(bs); rapi->LogOutput("%s", (pTempBuffer) ? pTempBuffer : "Unknown"); } break; case skVR_AttributeTag: { unsigned short tagValues[2]; elem.ReadArrayElementValues<unsigned short>(tagValues, 2, bs, transferSyntax); rapi->LogOutput("(%04x,%04x)", tagValues[0], tagValues[1]); } break; case skVR_Float: case skVR_Double: { const double value = elem.ReadSingleElementValue<double>(bs, transferSyntax); rapi->LogOutput("%f", value); } break; case skVR_SignedLong: case skVR_SignedShort: { const int value = elem.ReadSingleElementValue<int>(bs, transferSyntax); rapi->LogOutput("%i", value); } break; case skVR_UnsignedLong: case skVR_UnsignedShort: { const unsigned int value = elem.ReadSingleElementValue<unsigned int>(bs, transferSyntax); rapi->LogOutput("%i", value); } break; case skVR_OtherByte: case skVR_OtherFloat: case skVR_OtherWord: case skVR_Sequence: rapi->LogOutput("..."); break; default: rapi->LogOutput("N/A"); break; } rapi->LogOutput("\n"); } } switch (elem.GetElementId()) { case skPixelDataElement: { unsigned char *pTempBuffer = NULL; int decompressedSize = 0; currentImage.mpData = (unsigned char *)bs.GetBuffer() + elem.GetValueOfs(); currentImage.mDataLength = elem.GetValueLength(); currentImage.mDataVR = elem.GetVR(); bool convertVR = false; switch (transferSyntax) { case kDTS_RLE: pTempBuffer = decompress_rle(&decompressedSize, currentImage.mpData, bs.GetSize() - elem.GetValueOfs(), currentImage.mBitsPerChannel, rapi); convertVR = true; break; case kDTS_JpegBaseline: case kDTS_JpegExtended: case kDTS_JpegSpectral: case kDTS_JpegProgressive: case kDTS_JpegLossless: case kDTS_JpegLosslessFirstOrder: case kDTS_JpegLSLossless: case kDTS_JpegLSNearLossless: case kDTS_Jpeg2000LosslessReversible: case kDTS_Jpeg2000LosslessOrLossy: case kDTS_Jpeg2000Part2LosslessReversible: case kDTS_Jpeg2000Part2LosslessOrLossy: pTempBuffer = decompress_jpeg(&decompressedSize, currentImage.mpData, currentImage, transferSyntax, bs.GetSize() - elem.GetValueOfs(), rapi); convertVR = true; break; default: break; } if (pTempBuffer && decompressedSize > 0) { //use the decompressed result if available currentImage.mpData = pTempBuffer; currentImage.mDataLength = decompressedSize; } if (convertVR && currentImage.mBitsPerChannel > 0) { //choose an appropriate VR based on the bits per channel, assuming our provided //VR is meaningless in the current transfer syntax. //OtherWord is picked for anything less than 32 that isn't exactly 8, as that //path provides for n-bit reads. if (currentImage.mBitsPerChannel >= 32) { currentImage.mDataVR = skVR_OtherFloat; } else if (currentImage.mBitsPerChannel != 8) { currentImage.mDataVR = skVR_OtherWord; } else { currentImage.mDataVR = skVR_OtherByte; } } //if we don't have a photometric interpretation (some files totally lack it), try to figure one out. if (currentImage.mPhotometricInterpretation == kDPI_Unknown && currentImage.mBitsPerChannel > 0 && currentImage.mDataLength > 0) { const int sliceCount = std::max<int>(currentImage.mSliceCount, 1); const int pixelCount = currentImage.mColumnCount * currentImage.mRowCount * sliceCount; const int pixelBitsForOneChannel = pixelCount * currentImage.mBitsPerChannel; const int channelCount = (currentImage.mDataLength * 8) / pixelBitsForOneChannel; switch (channelCount) { case 1: currentImage.mPhotometricInterpretation = kDPI_Monochrome2; break; case 3: currentImage.mPhotometricInterpretation = kDPI_RGB; break; //still haven't found any argb in the wild to implement it. /* case 4: currentImage.mPhotometricInterpretation = kDPI_ARGB; break; */ default: //sorry, no idea break; } } //try converting the data, and reset the local image structure. //not sure if we're always guaranteed to have the data preceded by all image properties, //but this is the case in all of the data i've seen thus far. convert_dicom_image(noeTex, currentImage, transferSyntax, rapi); currentImage = SDicomImage(); if (pTempBuffer) { rapi->Noesis_UnpooledFree(pTempBuffer); } } break; case skPixelDataPlanarConfigurationElement: currentImage.mPlanarConfiguration = elem.ReadSingleElementValue<int>(bs, transferSyntax); break; case skPixelDataSamplesPerPixelElement: currentImage.mSamplesPerPixel = elem.ReadSingleElementValue<int>(bs, transferSyntax); break; case skPixelDataBitsAllocatedElement: currentImage.mBitsPerChannel = elem.ReadSingleElementValue<int>(bs, transferSyntax); break; case skPixelDataBitsStoredElement: currentImage.mBitsStoredPerChannel = elem.ReadSingleElementValue<int>(bs, transferSyntax); break; case skPixelDataHighBitElement: currentImage.mHighBit = elem.ReadSingleElementValue<int>(bs, transferSyntax); break; case skPixelDataPhotometricInterpretationElement: { elem.PrepDataStreamForValueRead(bs, transferSyntax); const char *pTempBuffer = elem.TemporaryStringForElementValue(bs); currentImage.mPhotometricInterpretation = photometric_interpretation_for_string(pTempBuffer); } break; case skPixelDataRowsElement: currentImage.mRowCount = elem.ReadSingleElementValue<int>(bs, transferSyntax); break; case skPixelDataColumnsElement: currentImage.mColumnCount = elem.ReadSingleElementValue<int>(bs, transferSyntax); break; case skPixelDataFrameCountElement: currentImage.mSliceCount = elem.ReadSingleElementValue<int>(bs, transferSyntax); break; case skPixelDataRescaleIntercept: currentImage.mRescaleIntercept = elem.ReadSingleElementValue<float>(bs, transferSyntax); break; case skPixelDataRescaleSlope: currentImage.mRescaleSlope = elem.ReadSingleElementValue<float>(bs, transferSyntax); break; case skRedPaletteLookupTableDescElement: set_palette_descriptor_from_element(¤tImage.mPaletteRgba[0], bs, transferSyntax, elem); break; case skGreenPaletteLookupTableDescElement: set_palette_descriptor_from_element(¤tImage.mPaletteRgba[1], bs, transferSyntax, elem); break; case skBluePaletteLookupTableDescElement: set_palette_descriptor_from_element(¤tImage.mPaletteRgba[2], bs, transferSyntax, elem); break; case skAlphaPaletteLookupTableDescElement: set_palette_descriptor_from_element(¤tImage.mPaletteRgba[3], bs, transferSyntax, elem); break; case skRedPaletteLookupTableDataElement: set_palette_data_from_element(¤tImage.mPaletteRgba[0], bs, elem); break; case skGreenPaletteLookupTableDataElement: set_palette_data_from_element(¤tImage.mPaletteRgba[1], bs, elem); break; case skBluePaletteLookupTableDataElement: set_palette_data_from_element(¤tImage.mPaletteRgba[2], bs, elem); break; case skAlphaPaletteLookupTableDataElement: set_palette_data_from_element(¤tImage.mPaletteRgba[3], bs, elem); break; default: break; } ++dataElementCount; elem.SetStreamToNextDataElement(bs); } if (logDataElements > 0) { rapi->LogOutput("Data Element count (may include sequence delimiters): %i\n", dataElementCount); } return true; } static void flush_group_stream(const unsigned short groupNum, const EDicomTransferSyntax outTS, RichBitStreamEx &outStream, RichBitStreamEx &groupStream) { //write the group length element CDicomDataElement(groupNum, 0x0000, skVR_UnsignedLong).WriteTypeToStream<unsigned int>( outStream, outTS, groupStream.GetOffset()); //then copy over all of the elements outStream.WriteBytes(groupStream.GetBuffer(), groupStream.GetOffset()); groupStream.SetOffset(0); } bool Image_MRI_SaveDicom(char *fileName, noesisTex_t *textures, int numTex, noeRAPI_t *rapi) { //potential optimization - lots of memory stream juggling, could be refactored RichBitStreamEx outStream; RichBitStreamEx groupStream; const int width = textures->w; const int height = textures->h; //only explicit little is supported at the moment const EDicomTransferSyntax outTS = kDTS_ExplicitVRLittleEndian; rapi->LogOutput("Exporting DICOM with transfer syntax %s...\n", skTransferSyntaxStrings[outTS]); for (int preambleIndex = 0; preambleIndex < skHeaderOffset; ++preambleIndex) { outStream.WriteByte(0); } outStream.WriteString("DICM"); unsigned char metaVerData[2] = { 0x00, 0x01 }; //file meta version CDicomDataElement(skGroup2, 0x0001, skVR_OtherByte).WriteToStream( groupStream, outTS, metaVerData, sizeof(metaVerData)); //media storage SOP - raw data CDicomDataElement(skGroup2, 0x0002, skVR_UID).WriteStringToStream( groupStream, outTS, "1.2.840.10008.5.1.4.1.1.6"); //media storage SOP instance UID CDicomDataElement(skGroup2, 0x0003, skVR_UID).WriteStringToStream( groupStream, outTS, "999.999.2.20150101.112000.2.666"); //transfer syntax UID CDicomDataElement(skGroup2, 0x0010, skVR_UID).WriteStringToStream( groupStream, outTS, skTransferSyntaxStrings[outTS]); //implementation class UID CDicomDataElement(skGroup2, 0x0012, skVR_UID).WriteStringToStream( groupStream, outTS, "999.999"); //end of group 2 flush_group_stream(skGroup2, outTS, outStream, groupStream); //image type CDicomDataElement(0x0008, 0x0008, skVR_CodeString).WriteStringToStream( groupStream, outTS, "ORIGINAL\\PRIMARY"); //SOP class UID CDicomDataElement(0x0008, 0x0016, skVR_UID).WriteStringToStream( groupStream, outTS, "1.2.840.10008.5.1.4.1.1.6"); //SOP instance UID CDicomDataElement(0x0008, 0x0018, skVR_UID).WriteStringToStream( groupStream, outTS, "999.999.2.20150101.112000.2.666"); //study date CDicomDataElement(0x0008, 0x0020, skVR_Date).WriteStringToStream( groupStream, outTS, "2015.01.01"); //content date CDicomDataElement(0x0008, 0x0023, skVR_Date).WriteStringToStream( groupStream, outTS, "2015.01.01"); //study time CDicomDataElement(0x0008, 0x0030, skVR_Time).WriteStringToStream( groupStream, outTS, "10:00:00"); //modality CDicomDataElement(0x0008, 0x0060, skVR_CodeString).WriteStringToStream( groupStream, outTS, "US"); //manufacturer CDicomDataElement(0x0008, 0x0070, skVR_LongString).WriteStringToStream( groupStream, outTS, "Noesis"); //referring physician's name CDicomDataElement(0x0008, 0x0090, skVR_PersonName).WriteStringToStream( groupStream, outTS, "Whitehouse, Dick"); //study description CDicomDataElement(0x0008, 0x1030, skVR_LongString).WriteStringToStream( groupStream, outTS, "Noesis export"); //series description CDicomDataElement(0x0008, 0x103E, skVR_LongString).WriteStringToStream( groupStream, outTS, "Very important data"); //stage name CDicomDataElement(0x0008, 0x2120, skVR_LongString).WriteStringToStream( groupStream, outTS, "Happy lovely time"); //stage number CDicomDataElement(0x0008, 0x2122, skVR_IntString).WriteStringToStream( groupStream, outTS, "1"); //number of stages CDicomDataElement(0x0008, 0x2124, skVR_IntString).WriteStringToStream( groupStream, outTS, "1"); //view number CDicomDataElement(0x0008, 0x2128, skVR_IntString).WriteStringToStream( groupStream, outTS, "1"); //number of views in stage CDicomDataElement(0x0008, 0x212A, skVR_IntString).WriteStringToStream( groupStream, outTS, "1"); //end of group 8 flush_group_stream(0x0008, outTS, outStream, groupStream); //patient full name CDicomDataElement(0x0010, 0x0010, skVR_PersonName).WriteStringToStream( groupStream, outTS, "Wong, Johnny"); //end of group 0x0010 flush_group_stream(0x0010, outTS, outStream, groupStream); //protocol name CDicomDataElement(0x0018, 0x1030, skVR_LongString).WriteStringToStream( groupStream, outTS, "Noesis Export"); //frame time if (numTex > 1) { CDicomDataElement(skFrameTimeElement, skVR_DecimalString).WriteStringToStream( groupStream, outTS, "33.333"); } //end of group 0x0018 flush_group_stream(0x0018, outTS, outStream, groupStream); //study instance UID CDicomDataElement(0x0020, 0x000D, skVR_UID).WriteStringToStream( groupStream, outTS, "999.999.2.20150101.112000"); //series instance UID CDicomDataElement(0x0020, 0x000E, skVR_UID).WriteStringToStream( groupStream, outTS, "999.999.2.20150101.112000.2"); //series number CDicomDataElement(0x0020, 0x0011, skVR_IntString).WriteStringToStream( groupStream, outTS, "2"); //instance number CDicomDataElement(0x0020, 0x0013, skVR_IntString).WriteStringToStream( groupStream, outTS, "666"); //end of group 0x0020 flush_group_stream(0x0020, outTS, outStream, groupStream); //samples per pixel CDicomDataElement(skPixelDataSamplesPerPixelElement, skVR_UnsignedShort).WriteTypeToStream<unsigned short>( groupStream, outTS, 3); //photometric interpretation CDicomDataElement(skPixelDataPhotometricInterpretationElement, skVR_CodeString).WriteStringToStream( groupStream, outTS, skPhotometricInterpretationStrings[kDPI_RGB]); //planar configuration CDicomDataElement(skPixelDataPlanarConfigurationElement, skVR_UnsignedShort).WriteTypeToStream<unsigned short>( groupStream, outTS, 0); //frame count and frame time increment pointer if (numTex > 1) { char tempString[256]; sprintf_s(tempString, "%i", numTex); CDicomDataElement(skPixelDataFrameCountElement, skVR_IntString).WriteStringToStream( groupStream, outTS, tempString); CDicomDataElement(skPixelDataFrameIncrementPointerElement, skVR_AttributeTag).WriteTypeToStream<unsigned int>( groupStream, outTS, skFrameTimeElement); } //row count CDicomDataElement(skPixelDataRowsElement, skVR_UnsignedShort).WriteTypeToStream<unsigned short>( groupStream, outTS, height); //column count CDicomDataElement(skPixelDataColumnsElement, skVR_UnsignedShort).WriteTypeToStream<unsigned short>( groupStream, outTS, width); //aspect ratio CDicomDataElement(skPixelDataAspectRatioElement, skVR_IntString).WriteStringToStream( groupStream, outTS, "4\\3"); //bits allocated CDicomDataElement(skPixelDataBitsAllocatedElement, skVR_UnsignedShort).WriteTypeToStream<unsigned short>( groupStream, outTS, 8); //bits stored CDicomDataElement(skPixelDataBitsStoredElement, skVR_UnsignedShort).WriteTypeToStream<unsigned short>( groupStream, outTS, 8); //high bit CDicomDataElement(skPixelDataHighBitElement, skVR_UnsignedShort).WriteTypeToStream<unsigned short>( groupStream, outTS, 7); //representation CDicomDataElement(skPixelDataRepresentationElement, skVR_UnsignedShort).WriteTypeToStream<unsigned short>( groupStream, outTS, 0); //end of group 0x0028 flush_group_stream(0x0028, outTS, outStream, groupStream); //build up an image buffer stream { RichBitStreamEx imageStream; for (int texIndex = 0; texIndex < numTex; ++texIndex) { noesisTex_t *pTexture = textures + texIndex; bool shouldFreeRgba = false; unsigned char *pRgba = rapi->Image_GetTexRGBA(pTexture, shouldFreeRgba); if (!pRgba) { rapi->LogOutput("WARNING: Could not retrieve RGBA data for texture %i.\n", texIndex); pRgba = (unsigned char *)rapi->Noesis_UnpooledAlloc(width * height * 4); memset(pRgba, 0, width * height * 4); shouldFreeRgba = true; } else { if (pTexture->w != width || pTexture->h != height) { //all image dimensions must match, resample it unsigned char *pResizedRgba = (unsigned char *)rapi->Noesis_UnpooledAlloc(width * height * 4); rapi->Noesis_ResampleImageBilinear(pRgba, pTexture->w, pTexture->h, pResizedRgba, width, height); if (shouldFreeRgba) { rapi->Noesis_UnpooledFree(pRgba); } //adopt the resized buffer after freeing the original buffer pRgba = pResizedRgba; shouldFreeRgba = true; } //write just the rgb values for (int pixelIndex = 0; pixelIndex < width * height; ++pixelIndex) { const unsigned char *pColor = pRgba + pixelIndex * 4; imageStream.WriteByte(pColor[0]); imageStream.WriteByte(pColor[1]); imageStream.WriteByte(pColor[2]); } } if (shouldFreeRgba) { rapi->Noesis_UnpooledFree(pRgba); } } //write the full image buffer to the group stream as pixel data CDicomDataElement(skPixelDataElement, skVR_OtherByte).WriteToStream( groupStream, outTS, imageStream.GetBuffer(), imageStream.GetOffset()); } //end of group 0x7FE0 flush_group_stream(0x7FE0, outTS, outStream, groupStream); const int finalSize = outStream.GetSize(); rapi->Noesis_WriteFile(fileName, outStream.GetBuffer(), finalSize); return true; } #define DICOM_DECL_OPTS(argRequired) \ dicomOpts_t *lopts = (dicomOpts_t *)store; \ assert(storeSize == sizeof(dicomOpts_t)); \ if (argRequired && !arg) \ { \ return false; \ } bool Model_MRI_FloatNormalizeHandler(const char *arg, unsigned char *store, int storeSize) { DICOM_DECL_OPTS(false); lopts->mFloatNormalize = true; return true; } bool Model_MRI_SBPHandler(const char *arg, unsigned char *store, int storeSize) { DICOM_DECL_OPTS(true); lopts->mApplySBP = true; sscanf(arg, "%f;%f;%f", &lopts->mSBPScale, &lopts->mSBPBias, &lopts->mSBPExponent); return true; } bool Model_MRI_DataNormalizeHandler(const char *arg, unsigned char *store, int storeSize) { DICOM_DECL_OPTS(false); lopts->mDataNormalize = true; return true; } bool Model_MRI_ApplySlopeHandler(const char *arg, unsigned char *store, int storeSize) { DICOM_DECL_OPTS(false); lopts->mApplySlope = true; return true; } bool Model_MRI_ElemLogHandler(const char *arg, unsigned char *store, int storeSize) { DICOM_DECL_OPTS(true); lopts->mLogDataElements = atoi(arg); return true; }
8313658 page hits since February 11, 2009.
Site design and contents (c) 2009 Rich Whitehouse. Except those contents which happen to be images or screenshots containing shit that is (c) someone/something else entirely. That shit isn't really mine. Fair use though! FAIR USE!
All works on this web site are the result of my own personal efforts, and are not in any way supported by any given company. You alone are responsible for any damages which you may incur as a result of this web site or files related to this web site.