[Insight-users] NiftiImageIO intensity rescaling bug
Tom Vercauteren
tom.vercauteren at m4x.org
Thu Mar 1 10:30:34 EST 2007
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