[Insight-users] Dicom Series Read Series Write: change in intensity values

[Matthias Schabel]

> I am getting kindda confused...but I'll try to answer your questions.
> Please let me know if I misunderstood you.
I'll try to clarify what I mean : I think that reading and then  
writing a DICOM file should do as little as necessary to either the  
DICOM data or the DICOM metadata. That is, if I read a file that has a  
specified rescaleSlope/Intercept and then write it out again, I should  
have a file with the same data and the same entries for rescaleSlope/ 
Intercept as initially.
> On Nov 27, 2007 2:09 AM, Matthias Schabel <mschabel at  
> ucair.med.utah.edu> wrote:
> >   That's correct MR also often -illegally- has a slope/intersept. So
> > does CT and I guess Secondary Capture. But we handle all of them
> > nicely when *reading* a DICOM file. The issue is that right now, we
> > are not able to store that image -just read- back to file.
What I don't understand is why you want to apply the rescaling terms  
during IO. Since rescaling is a peculiarity of DICOM files, users who  
read DICOM files need to know about their possible presence. Rescale  
slope and intercept are basically a prescription for how one should  
obtain a floating point value from the int16 values contained in the  
DICOM file, so they should only be used when converting an int16 DICOM  
data set to float/double - IO shouldn't be responsible for that at  
all. If you just ignore the presence of the rescales and pass them  
through on either read or write, you get the same data on write that  
you read in initially.
> >  To go into further details, PET for instance has a dynamic
> > slope/intersept. This make not only very efficient for this type of
> > data, but also completely necesseray since the scalar range is so  
> wide
> > and would prevent people f rom storing such images on 16bits as per
> > DICOM definition. Which means that ITK should have a mechanism
> > (callback?) which would allow GDCM to pass the new slope/intercept  
> for
> > each new 2D image when iterating over a 3D volume. So far I have not
> > look into it too much.
> >
> > But in summary: reading 'floating' point is handled, writing it back
> > is not unfortunately :(
Since DICOM has metadata for each slice, the rescaleSlope/Intercept  
can be different for each slice. When someone wants to get floating  
point values from DICOM data, what should happen, in my opinion, is :
1) DICOM data and metadata are read and stored.
2) User explicitly applies rescaling (maybe through  
itkDicomRescaleFilter<float/double>) to get floating point values.
> >  What I would expect is that either (in
> > order of my preference)
> >
> > 1) rescaling would only happen if explicitly invoked
>
> Per slice ? Per volume ? How would you know when to apply it in a
> batch environment ? the value loaded would not even have the correct
> unit for modality such as CT...
Per slice if metadata with rescale factors is provided for each slice,  
otherwise per volume, I guess. My point is that people need to know  
that the DICOM file doesn't contain Hounsfield Units for CT, just an  
integer discretization that approximates them. The precision of the  
conversion from the int16 values to HU will depend on the specific  
choice made for rescaling factors. If rescaleSlope is 0.1, then the HU  
increment will be 0.1, if it is 0.5 then the increment will be 0.5,  
but a larger range of HUs will be representable.
> > 2) rescaling would happen invertibly when DICOM images are  
> converted to
> > float and vice versa...
>
> I do not understand that sentence...
If you read a DICOM image, then convert to float by applying rescaling  
and, possibly, some other image processing operations on the data, and  
then want to write the resulting data back out to a new DICOM file,  
the float values can be converted back to int16 using the original  
rescaling factors (or new ones, if they are explicitly provided).
> > Since you can't correctly apply rescaling for integer conversions  
> without
> > loss of precision, it seems like a bad choice for the default  
> behavior. As
>
> Again I do not understand that sentence. Your Stored Pixel are 16bits,
> you apply an affine transform with floating point parameters (32bits)
> and load your real world value into 64bits double types... There is no
> way you can have any kind of loss of precision in this case, right ?
> Your input data has to be coming from some kind of acquisition device
> where it is dealing with descretized values, so the only loss I can
> think of is when this device is storing the discretized value as an
> approximation of an affine transform and a 16bits value.
Actually, single precision float can only represent integers up to  
2^24 exactly and double precision can only represent integers up to  
2^53 exactly. Thus, if you have a value of 32767 (largest  
representable in int16) and a rescale slope of anything larger than  
512 you won't be able to recover the exact initial integer (that is  
2^15*2^9 - (2^15-1)*2^9 will give you zero instead of 1). This may  
seem pedantic, but it is important to recognize the potential for  
errors here.
> > Those entries in the metadata seem to me to
> > be providing the additional information necessary to do integer- 
> >float and
> > float->integer conversions.
>
> We simply give up when the itk::Image is of type float...so there is
> no easy way to shoot yourself in the foot ...
> Removing the Rescale slope/Intercept is to work aroung the issue I
> described in a previous email, where input data is passed back to
> output.

I don't understand this comment - ITK's change in behavior to silently  
removing my metadata has shot me in the foot, while requiring explicit  
application of rescaling when processing DICOM files shouldn't shoot  
anyone in the foot if they know what they're doing. Basically, the  
data in the DICOM files for CT is not in HU and needs to be converted,  
the data in the DICOM files for PET is not in SUV and needs to be  
converted, etc... This should be done explicitly by the user, not in  
the ITK IO code, IMHO...
In any case, I can always regress back to ITK 3.0 if necessary, but I  
think it would be better to come up with a resolution to this issue.  
I'd love to hear other opinions on the matter.
Regards,
Matthias