Group Abstract Group Abstract

Message Boards Message Boards

Load DICOM with actual data values: Import/Export changes data

I'm very happy with the continued improvement of the Dicom import and export functionality and speed. However, there is one issue with the current implementation that makes it very unuseful if you do quantitative image analysis. It might be that I miss an option if not I think this should be fixed.

In the Mathematica documentation and examples, Dicom data is typically shown and imported as images which I understand for display purposes. But I consider the information in Dicom files as data, very well curated and standardized. For many (MRI) applications, the actual quantitative values of voxels stored in Dicom actually have meaning, values, and even units. For example in the image below each voxel value is actually a quantitative measure of T2 relaxation time in the heart, where the values are stored in milliseconds as voxel values as is also mentioned in the metadata.

enter image description here

enter image description here

To get from the stored values to quantitative values (WV, DV or FP) the fields from the header and the equations are well defined.

Header values:

  • SV = stored value of DICOM PIXEL DATA without scaling
  • WS = RealWorldValue slope (0040,9225) "RWVSlope"
  • WI = RealWorldValue intercept (0040,9224) "RWVIntercept"
  • RS = rescale slope (0028,1053) "RescaleSlope"
  • RI = rescale intercept (0028,1052) "RescaleIntercept"
  • SS = scale slope (2005,100E) "ScaleSlope"

Outputs:

  • WV = real world value
  • FP = precise value
  • DV = displayed value

Formulas:

  • WV = SV * WS + WI
  • DV = SV * RS + RI
  • FP = DV / (RS * SS)

So my first try was that I want to obtain the "RawData" to access the SV pixel data, which does not output anything.

enter image description here

Eventually, if I import this Dicom file into Mathematica I have to use a lot of tricks to get to the correct stored values. The Dicom images I use are stored as 12-bit Integers but are converted by Mathematica to a Numerical array with type Int16. Also, I have to specifically specify that I don't want any "DataTransformation" which by default rescales the data and actually changes some voxel values!!!!

In[1]:= << QMRITools`

{meta, data, bd} = 
  Import[file, {"dicom", {"MetaInformation", "Data", "BitDepth"}}, 
   "DataTransformation" -> None];
dataT = Import[file, {"dicom", {"Data"}}];

{dd = ToExpression[
   StringJoin @@ 
    StringCases[NumericArrayType[data], DigitCharacter]], bd}
{ss, rs, ri} = 
 meta /@ {"2005_100e", "RescaleSlope", "RescaleIntercept"}

{data, dataT} = Normal@{data, dataT};

svT = 2.^bd (dataT/(2.^dd));
pfT = (rs svT + ri)/(rs ss);

sv = 2.^bd (data/(2.^dd));
pf = (rs sv + ri)/(rs ss);

PlotData[pfT, pf]

Out[3]= {16, 12}

Out[4]= {1.99854, 0.500366, -2.}

For my current research, I need the PF values and to get them correctly I have to:

  1. Use "DataTransformation"->None, which is not really well documented what it actually does. But based on what I see actual values of the data are changed by clipping the histogram, which for default handling of medical data is never OK!! enter image description here
  2. As far as I am aware the only way to obtain the actual stored values of my Dicom data (the actual binary values stored in the file itself) I have to find the actual BiteDepth and the imported data type and rescale my data accordingly.

Below are the obtained PF valued data I need with and without DataTransformation. Although the image on the right might look less appealing with default range and scaling it is actually correct when scaled apropriately. enter image description here enter image description here

Am I missing a correct option for getting the actual data stored? If not this should definitely be changed. Dicom is the international standard to transmit, store, retrieve, print, process, and display medical imaging information. With the current implementation, it is impossible to Import and Export such files without actually changing the stored data.

Thanks, Martijn

POSTED BY: Martijn Froeling
5 Replies

I am coming upon this interesting thread years later. I want to comment on the sentence, "We went for this design because otherwise a vast majority of DICOM files would return an all-black image by default, which wouldn't be very useful." I think there's a hidden assumption about to whom the functionality is meant to be useful.

With automatic transformation and glossy images, Import[]'s DICOM handling is useful to non-professionals who want to quickly put an existing DICOM image in a format that is easy to qualitatively analyze. That may be useful to patients who are generating informed questions for their healthcare providers. It is also useful to people who are running demos of Wolfram's functionality to generate sales.

With automatic transformation and glossy images, Import[]'s DICOM handling is NOT useful to healthcare providers, because as the OP mentioned: "...based on what I see actual values of the data are changed by clipping the histogram, which for default handling of medical data is never OK!!" Conventions like this are very serious in medicine because doctors need images free of manipulations that could affect the interpretation of key features of the image. So anyone who wanted to build a medical image processing app in Mathematica that meets government standards would have to show they are explicitly working around this default behavior.

Without automatic transformation, Importing a DICOM indeed gets you mostly-black images initially, but it seems that is what's useful to researchers. The OP was a very polite way of saying that without the ability to extract the raw voxel values, Wolfram's exciting new Import functionality for DICOMs was effectively useless to researchers. And the solution was for researchers to use a then-undocumented "RawPixelData" tag. (There's a single line about it in this article https://reference.wolfram.com/language/ref/format/DICOM.html now... was that there 4 years ago?) But in research, reproducibility is paramount, and reproducibility is curtailed greatly by "black box" functionality in software.

I wonder if sometimes the community would be better served by more candor and clarity about which functionality is meant to be useful to which users. Right now the main documentation page for DICOM reads like it's geared toward people who have never heard of a DICOM image until that day, and it's very important to keep that language in place in order to continue facilitating learning. But the info a researcher needs to do even 1 micro-step of their work ("RawPixelData" --> "DICOM pixel data without scaling") is buried in the Import Elements section, and a person who doesn't read that far is making a rational choice since Wolfram documentation just doesn't always cover the relevant technical information. This post was a demonstration of someone having to figure out for themself what's happening "under the hood" in order to make the functionality useful for their application area.

A fix could possibly look like creating an "Application Areas" section on documentation pages. For the DICOM page (https://reference.wolfram.com/language/ref/format/DICOM.html) that section would be close to the top and it would say something right at the beginning like:

"Import[] by default scales pixel values by clipping the histogram to ensure that the first image generated from the data has visible features. Researchers who want to perform a quantitative analysis must Import the "RawPixelData" element or Import with "DataTransformation" -> None".

Doing this would potentially allow researchers to get to their real work faster, without compromising other people's learning.

Awesome! Very helpful information!

I was not aware of the DICOMTools will dig a bit into that.

And indeed the 2005 tag is a vendor-specific one, if i remember correctly all odd first numbers are generally vendor-specific.

Thanks again for the information!

POSTED BY: Martijn Froeling

Importing a list of specific tags in one call to Import is not possible right now, but we will try to add support for this in the next release. Currently, MetaInformation can only handle one subelement which must be a string. In the simplest case this will be a name of a DICOM tag, e.g. "PixelData", but it is possible to pass a simple pattern to import a bunch of related tags at once, for instance:

In[94]:= Import["ExampleData/head.dcm.gz", {"MetaInformation", "PatientName"}]

Out[94]= "Male 1958 Anonymous"

In[95]:= Import["ExampleData/head.dcm.gz", {"MetaInformation", "Patient*"}]

Out[95]= <|"PatientName" -> "Male 1958 Anonymous", 
 "PatientID" -> "GH-WISE:68736", "PatientSex" -> "M", 
 "PatientBirthName" -> "anonymous", 
 "PatientMotherBirthName" -> "anonymous", 
 "PatientTelephoneNumbers" -> "anonymous", 
 "PatientReligiousPreference" -> "anonymous", 
 "PatientComments" -> "anonymous", "PatientPosition" -> "HFS"|>

That being said, you can still read a list of tags without calling Import multiple times if you are willing to use DICOMTools paclet directly (Import uses this paclet under the hood). Paclet functions are not documented and we don't guarantee backwards compatibility for them, but they often offer more flexibility then Import interface:

 In[97]:= Needs["DICOMTools`"]

 In[98]:= file = First @ ExtractArchive @ FindFile @ "ExampleData/head.dcm.gz";

 In[99]:= DICOMTools`ReadMetadata[file, {"PatientName", "RescaleSlope", "RescaleIntercept"}]

 Out[99]= <|"PatientName" -> "Anonymous^Male 1958",  "RescaleSlope" -> "1.53968253968253", "RescaleIntercept" -> "0.0"|>

Two things to notice: paclet only works with full paths and does not handle compressed DICOMs, so we need an extra step to uncompress the sample file and secondly, values are imported as they are stored in the file, without post-processing or formatting (so for instance, numbers are imported as strings and dates will not be automatically converted to WolframLanguage DateObjects).

If you wish to format the returned data, another step is needed:

In[113]:= rawMeta = DICOMTools`ReadMetadata[file, {"PatientName", "RescaleSlope", "RescaleIntercept"}]

Out[113]= <|"PatientName" -> "Anonymous^Male 1958",  "RescaleSlope" -> "1.53968253968253", "RescaleIntercept" -> "0.0"|>

In[114]:= DICOMTools`ConvertMetadata[file][rawMeta]

Out[114]= <|"PatientName" -> "Male 1958 Anonymous", "RescaleSlope" -> 1.5396825, "RescaleIntercept" -> 0.|>

To answer your second question - Import only supports tag keywords, but the paclet can read tags using (group, elem) values. They can be specified in the form: {"GGGG", "EEEE"} (leading 0s can be omitted). For example:

In[105]:= DICOMTools`ReadMetadata[file, {{"0028", "1053"}, {"10", "40"}}]

Out[105]= <|"RescaleSlope" -> "1.53968253968253", "PatientSex" -> "M"|>

In the returned Association tag names are used as keys anyway, but if you prefer, you can always replace any tag name with its group-elem id using the following function:

In[107]:= DICOMTools`GetGroupElemNumbers["PixelData"]

Out[107]= {"7FE0", "10"}

In[108]:= DICOMTools`GetGroupElemNumbers["PatientAge"]

Out[108]= {"10", "1010"}

There is also a reverse function:

In[109]:= DICOMTools`GetKeyword[{"7FE0", "0010"}]

Out[109]= "PixelData"

In[110]:= DICOMTools`GetKeyword[{"10", "1010"}]

Out[110]= "PatientAge"

Notice that it only works with standard DICOM tags. I guess that 2005_100e that you use in your example is vendor-specific and so the paclet will not know its name:

In[111]:= DICOMTools`GetKeyword[{"2005", "100e"}]

Out[111]= "Unknown Tag & Data"
POSTED BY: Rafal Chojna

Awesome, that was exactly what I was looking for!!!

Thanks for the explanations, always nice to learn a few new things, I was not aware that with {"MetaInformation", "xxx"} you can only import specific tags.

Two additional questions if i may:

  1. Is loading multiple tags in one go possible based on the help of Import this should work right? enter image description here

These examples do not work or am I again missing something.

    Import["data\\dcm\\001_v1\\01901_T1_enhanced_(3x8mm)\\00001.dcm", \
    {"MetaInformation", "2005_100e", "RescaleSlope", "RescaleIntercept"}]
    Import["data\\dcm\\001_v1\\01901_T1_enhanced_(3x8mm)\\00001.dcm", \
    {"dicom", "MetaInformation", "2005_100e", "RescaleSlope", 
      "RescaleIntercept"}]
  1. For known Dicom tags only the Name is exposed, but not the orriginal tag, is this correct? For exmaple Pixel Data has tag (7FE0,0010)

    (*works*)
    Import["data\\dcm\\001_v1\\01901_T1_enhanced_(3x8mm)\\00001.dcm", \
    {"MetaInformation", "PixelData"}]
    (*does not work*)
    Import["data\\dcm\\001_v1\\01901_T1_enhanced_(3x8mm)\\00001.dcm", \
    {"MetaInformation", "7FE0_0010"}]
    

Thanks again for the support!

Best Martijn

PS Dicom support came a long way, I remember using MathVisionTools for Dicom support back in my university days

POSTED BY: Martijn Froeling

Hi Martijn,

I'm sorry to hear that you are having a hard time working with DICOM files in the Wolfram Language. As you correctly noticed, Import["f.dcm", "Image"] applies a number of transformations under the hood. We went for this design because otherwise a vast majority of DICOM files would return an all-black image by default, which wouldn't be very useful.

I admit that our documentation for DICOM could be improved to be more helpful for advanced users. The element you are looking for is called "RawPixelData" and for now remains undocumented.

In[41]:= rawData = Import["ExampleData/head.dcm.gz", "RawPixelData"]
Out[41]= NumericArray[Type: UnsignedInteger16 Dimensions: {512, 512}]

In[42]:= MinMax[Normal @ rawData]
Out[42]= {0, 1727}

It returns an array of 16-bit unsigned integers even when the actual bit depth of the values in the file is 12, but the values are not modified in any way.

When you Import the "Data" element, on the other hand, you get the fully processed data which corresponds to the pixel values of the "Image" element:

enter image description here

Now, with "DataTransformation" -> None you can import pixel data without any scaling or other transformations defined in the DICOM file, but it will still map the theoretical range of pixel values to the range of the NumericArray type. For instance, a 12-bit DICOM data will be returned in a NumericArray of type "UnsignedInteger16" (because we don't have a 12-bit wide in the system) and every value will be multiplied by 16 (which maps [0, 2^12) to [0, 2^16)):

In[78]:= dataNoTransform = Import["ExampleData/head.dcm.gz", "Data", "DataTransformation" -> None]
Out[78]= NumericArray[Type: UnsignedInteger16 Dimensions: {512, 512}]

In[79]:= MinMax[Normal @ dataNoTransform]
Out[79]= {0, 27632}

In[80]:= MinMax[Normal @ rawData]
Out[80]= {0, 1727}

In[81]:= Max[dataNoTransform] / Max[rawData]
Out[81]= 16

One last tip: every standard DICOM tag can be imported via the "MetaInformation" element, including pixel data:

In[85]:= pixData = First @ Import["ExampleData/head.dcm.gz", {"MetaInformation", "PixelData"}]
Out[85]= NumericArray[Type: UnsignedInteger16 Dimensions: {512, 512}]

In[86]:= pixData === rawData
Out[86]= True
POSTED BY: Rafal Chojna
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard