[Insight-users] Resampling an image to iso voxels

[Miller, James V (Research)]

The transformation supplied to the resample filter maps the physical 
coordinates (mm) of a pixel in the output image to the physical 
coordinates (mm) of a pixel in the input image. (You can look at pages
178-185 of the Software Guide for more details.)

In your example, it appears as though you want not only isotropic 
voxels but voxels that are 1mm on a side. (You could have isotropic
voxels that have any length on a side if you wanted. For instance, suppose
the
input spacing is 0.4x0.4x1.5, then I would have resampled the image to 
be 0.4x0.4x0.4 so only the resolution in z is changed. My resulting pixels 
are isotropic, just not unit volume.).

Since the transformation maps output physical coordinates to input
physical coordinates, I think you want the "scales" to be the input
spacing and not 1.0 / input spacing.

Suppose the input spacing is 1.5 and let's look at pixel 
index (1, 1) in the output image, this pixel is at is (1mm, 1mm) 
in physical coordinates since the output spacing is 1.  If the scales 
are set are you specified as 1/1.5 = .666, then the corresponding pixel 
in the input image is at physical coordinates (.666 mm, .666 mm).  To 
convert this to an index, we divide by the spacing or 1.5.  This gives us 
an index of (.444px, .444px).  This is what you saw.

If we set the scales to 1.5, then the corresponding pixel in the input
image is at physical coordinates (1.5mm, 1.5mm). To convert this to an
index, 
we divide by the spacing or 1.5. This gives us an index of (1px, 1px).
Which is probably what you want.

Now if we take an image that as spacing 1.5, use 1.5 as the scales on the 
transformation and set the output spacing to be 1.0. Then our output image 
will have the same number of pixels as our input image. All we have ended up
doing is changing the spacing, we haven't needed to interpolate any pixels.
If we set the scales to 1.0 and set the output spacing to be 1.0, then our 
output image must be 1.5 times as big and we will interpolate pixels 
as needed. I think this is more in line with your intentions.

(If our input spacing was 0.5, output spacing of 1.0, and scales of 1, 
then out output image will be half the size of the input.)

Jim


-----Original Message-----
From: Ralf o Floca [mailto:rfloca at stud . fh-heilbronn . de]
Sent: Thursday, December 04, 2003 3:55 AM
To: insight-users at itk . org
Subject: [Insight-users] Resampling an image to iso voxels



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