[Insight-users] Resampling an image to iso voxels

[Luis Ibanez]

Hi Ralf,

For converting your image into an isotropic image
you don't need to use the scaling transform, just use
the *Identity* transform.  Since ITK resampling happens
in physical coordinates, and transforms specify the
convertion from physical coordinates A to physical
coordinates B, you simply need to set a null transform.

The actual resampling is performed by specifying a new
*pixel spacing* for the output image as well as a new
image size.

You will find a new example on how to convert images
to Isotropic voxels under

   Insight/Examples/Filtering/
                ResampleVolumesToBeIsotropic.cxx


http://www . itk . org/cgi-bin/cvsweb . cgi/Insight/Examples/Filtering/ResampleVolumesToBeIsotropic . cxx?cvsroot=Insight&sortby=date

You actually use this example as a command line utility.

If you look at the code in this example, please note that
you must smooth an image along the dimensions that are going
to be sub-sampled. This is needed to avoid the superposition
of spectra in the Fourier space.


Please let us know if you have any problems,


Thanks


   Luis


---------------------
Ralf o Floca wrote:

> Hello
> 
> I have a problem with resampling in the following situation:
> I have a few 3d images (analyze or stacks of 2D dicom images) with a
> spacing in some cases unequal to 1.
> After reading the volumes, I want to convert them to an iso voxel
> spacing of 1. I have chosen a resample filter with a scale transform to
> manage this task, the code for this is something like that:
> 
> void ConvertImageToIsoSpacing(OutputImageType::Pointer& m_Image)
> {
>   itk::FixedArray< double, iDimension> scales;
>   const double* pdSpacing = m_Image->GetSpacing();
>   double dNewSpacing[3];
> 
>   itk::Size<iDimension> newSize;
>   const itk::Size<iDimension>& oldSize =
> m_Image->GetLargestPossibleRegion().GetSize();
> 
>   for (int iIndex=0; iIndex< 3; iIndex ++, pdSpacing++)
>   { //Calculate the scaling for the transform and the size of the output
>     dNewSpacing[iIndex] = 1.0;
>     if (*pdSpacing!=0) scales[iIndex] = 1/(*pdSpacing);
> 		
>     newSize[iIndex] = oldSize[iIndex]*(*pdSpacing);
>   };
> 
>   OutputImageType::Pointer m_ImageTemp = OutputImageType::New();
> 
>   typedef itk::ResampleImageFilter< OutputImageType, OutputImageType >
> ResampleType;
>   typedef itk::ScaleTransform< double, iDimension > ScaleType;
>   
>   ResampleType::Pointer resampler = ResampleType::New();
> 		
>   ScaleType::Pointer scaler = ScaleType::New();
>   scaler->SetScale(scales);
> 	  
>   resampler->SetTransform(scaler);
>   resampler->SetSize(newSize);
>   resampler->SetOutputOrigin(m_Image->GetOrigin());
>   resampler->SetOutputSpacing(dNewSpacing);
>   resampler->SetDefaultPixelValue( 0 );
>   resampler->SetInput(m_Image);
> 
>   m_ImageTemp = resampler->GetOutput();
>   resampler->Update();
>   m_Image = m_ImageTemp;
> };
> 
> For dicom-stacks it seems to work fine, but when using it with an image
> saved as an analyze the resampling produces bad results for spacings
> unequal to 1, because the output images seemed to be double scaled
> (e.g.: when the spacing is 1.5, the image is scaled by 0.444.. -and not
> 0.666..-, but the new image size is correct, so overlapping parts are
> cut away). I was able to locate the reason for the "double use" of the
> scaling, but I don't know where the misuse is in my code.
> 
> 1. 1st time the scaling is used is in the resample filter
> (itkResampleImageFilter.txx, Line 220):
> 	// Compute corresponding input pixel position
> 	inputPoint = m_Transform->TransformPoint(outputPoint);
> There the outputPoint is scaled the first time.
> 
> 2. Then by the interpolation in line 227, which uses
> itkImage::TransformPhysicalPointToContinousIndex() (itkImage.h, line
> 259), where the inputPoint origin difference is divided by the input
> image spacing. At the end the output index has been divided twice by
> spacing to get the input index.
> 
> If it not an error in the design, I would be glad if someone can show me
> my misuse of itk.
> 
> The reason, why the dicom stacks are converted without error in most
> cases, leads me to my next question/problem. I extrapolate the z spacing
> myself from the dicom stack and actualize the image, loaded with the
> help of the itkImageSeriesReader. To do this, I use following piece of
> code:
> 
> 	const double* pdSpacing = m_Image->GetSpacing();
> 	double dNewSpacing[3] = {1};
> 
> 	for (int iIndex=0; iIndex< 2; iIndex ++, pdSpacing++)
> 	{
> 	  dNewSpacing[iIndex] = *pdSpacing;
> 	};
> 
> 	dNewSpacing[2] = dSliceThickness 
> 	//dSliceThickness is the z spacing extrapolated from the dicom
> slices;
> 
> 	m_Image->SetSpacing(dNewSpacing);
> 
> The spacing is now actualized and every thing seems to be fine. But when
> I convert the image (with the function listed above) and reach the
> BeforeThreadedGenerateData function of the ResampleImageFilter (line
> 174) where the input image of the interpolator is to be set, the spacing
> of the input image is 1.
> So there is no error when converting the image, because the input point
> is divided by 1 when converting it to index, but I don't know, what I'm
> doing wrong to lose the spacing information.
> 
> Thank you, for working through this long mail. Hope, someone can give me
> a or some clues.
> 
> Ralf o Floca
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk . org
> http://www . itk . org/mailman/listinfo/insight-users
>