[Insight-users] Resampling an image to iso voxels

[Ralf o Floca]

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