[Insight-users] NiftiImageIO intensity rescaling bug
Hans Johnson
hans-johnson at uiowa.edu
Thu Mar 8 09:13:32 EST 2007
Tom,
Thank you for all your work in identifying both the problem and a potential
solution. The issues you bring up here, and your thorough investigation of
the problem are certainly very helpful.
ITK is in the middle of preparing for a release, and so these types of
changes will not be accepted right now. It will be 10-20 days before this
issue can really be investigated fully.
A gut reaction so far is that the rescaling proposal is reasonable and
logical for making ITK compliant to the stated NIFTI format intention. The
ConvertPixelBuffer seems to be an over-specialization to address your issue,
and seems like it is not generalizable in all cases. Would it be possible
for the you to read in the nifti segmentation images as floating point, and
then do the ConvertGrayToGray behavior explicitly outside of the file reader
(thus pushing you needed specialization out of the general purpose IO
mechanism)?
If you have not done so already, please submit a bug report so that when
this problem is addressed there will be a permanent log of what the
resolution was.
Hans
On 3/8/07 5:31 AM, "Tom Vercauteren" <tom.vercauteren at m4x.org> wrote:
> Hi Hans,
>
> After some more thought, I think that my previous solution needed more
> work. I found that I needed to modify both the NiftiImageIO and the
> ConvertPixelBuffer classes.
>
> This seems a little involved by I didn't find any better solution.
> Please tell me I am fooling myself.
>
>
> 1) What I did to modify NiftiImageIO
> ==========================
> Here is what I understand it from the NIFTI documentation:
> http://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_page
> s/scl_slopeinter.html
>
> We need to rescale the data if scl_slope!=0 or (scl_slope!=1 and
> scl_inter!=0). If we need to rescale the data, then the pixel type
> should be promoted to a floating point data type. It should be "float"
> if the written pixel type is not a floating point data type. Otherwise
> it can be left unchanged.
>
> In case we need to rescale, we should thus:
> - determine the pixel type based on the written pixel type
> - read the data into a written pixel type buffer
> - rescale the written pixel to the final pixel type buffer
>
> I have attached my modified NIFTIImageIO class. It is not completely
> finished yet as I always use floats if a rescaling is needed but it
> gives an idea of how this can be done. I hope this helps.
>
>
> 1) What I did to modify ConvertPixelBuffer
> ==============================
> But If I only do that I am still stuck with the roundoff error problem
> I mentioned. The ImageIO class now provides floating point data. If I
> ask the ImageFileReader to read a integer image, it is up to the
> ConvertGrayToGray function in itkConvertPixelBuffer.txx to do the
> conversion job. Here is the problem:
> - it uses a static_cast to convert the pixel type provided by the
> ImageIO to the pixel type asked by the ImageFileReader
> - If the ImageIO pixel type is a floating point type and the
> ImageFileReader pixel type is an integer type, then this leads to
> round-off errors. These errors can be really problematic when one
> reads segmentation images for example
> - the solution I found was to use a rounding operation before the
> cast. See below for the ConvertGrayToGray function I use.
> template < typename InputPixelType,
> typename OutputPixelType,
> class OutputConvertTraits
>>
> void
> ConvertPixelBuffer<InputPixelType, OutputPixelType, OutputConvertTraits>
> ::ConvertGrayToGray(InputPixelType* inputData,
> OutputPixelType* outputData , int size)
> {
> InputPixelType* endInput = inputData + size;
>
> if ( NumericTraits<OutputComponentType>::is_integer )
> {
> while(inputData != endInput)
> {
> OutputConvertTraits
> ::SetNthComponent(0, *outputData++,
> static_cast<OutputComponentType>
> (vnl_math_rnd(static_cast<double>(*inputData))));
> inputData++;
> }
> }
> else
> {
> while(inputData != endInput)
> {
> OutputConvertTraits
> ::SetNthComponent(0, *outputData++,
> static_cast<OutputComponentType>
> (*inputData));
> inputData++;
> }
> }
> }
>
>
> Best,
> Tom
>
>
>
>
>
> On 3/1/07, Tom Vercauteren <tom.vercauteren at m4x.org> wrote:
>> Hans,
>>
>> I found another problem with the rescaling function. It uses a cast to
>> convert the the double to the chosen PixelType. If the PixelType is an
>> integer this leads to roundoff errors.
>>
>> To make things clear, I will briefly mention what I need to do. I am
>> trying to correctly read the images from BrainWeb
>> (http://www.bic.mni.mcgill.ca/brainweb/anatomic_normal_20.html).
>>
>> The images are stored in MINC format which is not (yet) natively
>> supported by ITK. I have therefore installed the MINC tools from
>> http://packages.bic.mni.mcgill.ca/tgz/ (instructions at
>> http://www.bic.mni.mcgill.ca/software/minc/ )
>>
>> I dowloaded a discrete segmentation in MINC format here:
>> http://www.bic.mni.mcgill.ca/cgi/brainweb1?alias=subject04_crisp&download=1
>>
>> I convert this MINC dataset into NIFTI with the mnc2nii tool:
>> mnc2nii -nii -short subject04_crisp_v.mnc.gz
>>
>> I then read the NIFTI image in ITK using a PixelType=short (I need to
>> use integers here because it is a segmentation and floating points
>> classes do not really make sense).
>>
>> When I first looked at the data I realized that some classes were
>> empty. Because of the rescaling, for the pixels where I should have a
>> I(pix)=3, I have something like 2.99 which results in a wrong tissue
>> classification I(pix)=static_cast<short>(2.99)=2.
>>
>> Here is what I did to correct this. In itkNiftiImageIO.cxx, I replaced
>> the RescaleFunction by 2 distinct functions:
>>
>> // Internal function to rescale pixel according to Rescale Slope/Intercept
>> template<class TBuffer>
>> void RescaleFunctionCast(TBuffer* buffer, double slope, double
>> intercept, size_t size)
>> {
>> for(unsigned int i=0; i<size; i++)
>> {
>> double tmp = static_cast<double>(buffer[i]) * slope;
>> tmp += intercept;
>> buffer[i] = static_cast<TBuffer>(tmp);
>> }
>> }
>>
>> template<class TBuffer>
>> void RescaleFunctionRound(TBuffer* buffer, double slope, double
>> intercept, size_t size)
>> {
>> for(unsigned int i=0; i<size; i++)
>> {
>> double tmp = static_cast<double>(buffer[i]) * slope;
>> tmp += intercept;
>> buffer[i] = static_cast<TBuffer>(vnl_math_rnd(tmp));
>> }
>> }
>>
>> And in the switch I call RescaleFunctionRound for integer types and
>> RescaleFunctionCast for floating point types:
>>
>> switch(m_ComponentType)
>> {
>> case CHAR:
>> RescaleFunctionRound(static_cast<char *>(buffer),
>> m_RescaleSlope,
>> m_RescaleIntercept,numElts);
>> break;
>> case UCHAR:
>> RescaleFunctionRound(static_cast<unsigned char *>(buffer),
>> m_RescaleSlope,
>> m_RescaleIntercept,numElts);
>> break;
>> case SHORT:
>> RescaleFunctionRound(static_cast<short *>(buffer),
>> m_RescaleSlope,
>> m_RescaleIntercept,numElts);
>> break;
>> case USHORT:
>> RescaleFunctionRound(static_cast<unsigned short *>(buffer),
>> m_RescaleSlope,
>> m_RescaleIntercept,numElts);
>> break;
>> case INT:
>> RescaleFunctionRound(static_cast<int *>(buffer),
>> m_RescaleSlope,
>> m_RescaleIntercept,numElts);
>> break;
>> case UINT:
>> RescaleFunctionRound(static_cast<unsigned int *>(buffer),
>> m_RescaleSlope,
>> m_RescaleIntercept,numElts);
>> break;
>> case LONG:
>> RescaleFunctionRound(static_cast<long *>(buffer),
>> m_RescaleSlope,
>> m_RescaleIntercept,numElts);
>> break;
>> case ULONG:
>> RescaleFunctionRound(static_cast<unsigned long *>(buffer),
>> m_RescaleSlope,
>> m_RescaleIntercept,numElts);
>> break;
>> case FLOAT:
>> RescaleFunctionCast(static_cast<float *>(buffer),
>> m_RescaleSlope,
>> m_RescaleIntercept,numElts);
>> break;
>> case DOUBLE:
>> RescaleFunctionCast(static_cast<double *>(buffer),
>> m_RescaleSlope,
>> m_RescaleIntercept,numElts);
>> break;
>> default:
>> if(this->GetPixelType() == SCALAR)
>> {
>> ExceptionObject exception(__FILE__, __LINE__);
>> exception.SetDescription("Datatype not supported");
>> throw exception;
>> }
>> }
>>
>>
>> I have added this to the bug I have filed previously. Hope this helps.
>>
>> Best,
>> Tom
>>
>> On 3/1/07, Hans J. Johnson <hans-johnson at uiowa.edu> wrote:
>>> Tom,
>>>
>>> If you can build a test case that fails with the current code, that would be
>>> the best way to determine what the correct limits are. I will not be able
>>> to write a test case in the near future, so your assistance would be greatly
>>> appreciated.
>>>
>>> Thanks,
>>> Hans
>>>
>>>
>>>
>>> On 2/26/07 3:46 AM, "Tom Vercauteren" <tom.vercauteren at m4x.org> wrote:
>>>
>>>> Hi Hans and Luis,
>>>>
>>>> Thanks for the fast patch!
>>>> Just one little thing. You chose to use
>>>> vcl_numeric_limits<double>::epsilon() to check the equality. That may
>>>> be a bit harsh. I don't know exactly how these numbers are stored on
>>>> the disk (float or double) but I guess that a more permissive check
>>>> using something like 4.0*vcl_numeric_limits<float>::epsilon() (at
>>>> least a factor 2 to account for the minus operation) would be
>>>> sufficient. What do you think?
>>>>
>>>> Best,
>>>> Tom
>>>>
>>>> On 2/23/07, Hans Johnson <hans-johnson at uiowa.edu> wrote:
>>>>> Tom,
>>>>>
>>>>> I agree with you. I will look more closely at this later today, and
>>>>> should
>>>>> have a cvs commit in before end of business today.
>>>>>
>>>>> Thank you for the detailed report and suggestion for a fix.
>>>>>
>>>>> Hans
>>>>>
>>>>>
>>>>> On 2/23/07 8:24 AM, "Tom Vercauteren" <tom.vercauteren at gmail.com> wrote:
>>>>>
>>>>>> Hi,
>>>>>>
>>>>>> I think there is a bug in the NiftiImageIO.
>>>>>>
>>>>>> I have a Nifti data set that has the following tags: scl_inter=0 and
>>>>>> scl_slope=0.0431. The data should be between 0 and 11 but when I read
>>>>>> it with ITK, I get values beween 0 and 255. The intensity rescaling is
>>>>>> thus not performed.
>>>>>>
>>>>>> I have been looking at the NiftiImageIO code and saw that the
>>>>>> intensity rescaling is only done if the following condition is met:
>>>>>> (m_RescaleSlope > 1 || m_RescaleIntercept != 0)
>>>>>>
>>>>>> I would have used something like:
>>>>>> ( fabs(m_RescaleSlope-1.0)>epsilon or fabs(m_RescaleIntercept)>epsilon )
>>>>>>
>>>>>> I have filed a bug report for this:
>>>>>> http://www.itk.org/Bug/bug.php?op=show&bugid=4489
>>>>>>
>>>>>> Best regards,
>>>>>> Tom Vercauteren
>>>>>> _______________________________________________
>>>>>> Insight-users mailing list
>>>>>> Insight-users at itk.org
>>>>>> http://www.itk.org/mailman/listinfo/insight-users
>>>>>
>>>>>
>>>
>>>
>>
More information about the Insight-users
mailing list