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, &params) != 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, &params) != 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(&currentImage.mPaletteRgba[0], bs, transferSyntax, elem);
			break;
		case skGreenPaletteLookupTableDescElement:
			set_palette_descriptor_from_element(&currentImage.mPaletteRgba[1], bs, transferSyntax, elem);
			break;
		case skBluePaletteLookupTableDescElement:
			set_palette_descriptor_from_element(&currentImage.mPaletteRgba[2], bs, transferSyntax, elem);
			break;
		case skAlphaPaletteLookupTableDescElement:
			set_palette_descriptor_from_element(&currentImage.mPaletteRgba[3], bs, transferSyntax, elem);
			break;
		case skRedPaletteLookupTableDataElement:
			set_palette_data_from_element(&currentImage.mPaletteRgba[0], bs, elem);
			break;
		case skGreenPaletteLookupTableDataElement:
			set_palette_data_from_element(&currentImage.mPaletteRgba[1], bs, elem);
			break;
		case skBluePaletteLookupTableDataElement:
			set_palette_data_from_element(&currentImage.mPaletteRgba[2], bs, elem);
			break;
		case skAlphaPaletteLookupTableDataElement:
			set_palette_data_from_element(&currentImage.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;
}


































































2823390 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.