[Insight-users] Problem with itk::ResampleImagefilter

[Luis Ibanez]

Hi Vincent,

Please take a look  at the diagram in the following Wiki page:

http://www.itk.org/Wiki/Proposals:Orientation#DICOM_LPS_Orientation_-_Full_Body_Graphical_Example

This diagram illustrates the default LPS coordinate system
used in ITK (following the DICOM standard).

When an image is loaded in ITK, and its "Direction Matrix"
is an Identity, then we expect the data set to match an LPS
orientation, and that will correspond to an "Axial" acquisition.

ITK images (since ITK 3.14) manage their internal orientation
correctly (that is, an itk::Image behaves just the same as an
itk::OrientedImage).

The final effect is that, if you want to extract an "Axial" slice
from an ITK image, you simply need to set the Output Direction
to be an Identity matrix

      *regardless of what the input image orientation is*

(which is a nice thing).


Now, if you need a Coronal Slice, then you should set the
Output Direction to be the matrix:

1   0   0
0   0  -1
0   1   0

and if you want a Sagittal slice at the output, then you
need to set the Output Direction to be the matrix:

0   1   0
0   0  -1
-1  0   0


This should take care of the orientation.

The second piece that you should take into account
is the origin of the Slice.

This is most likely the source of the problem that you
are encountering.

You should look at the figure
http://www.itk.org/Wiki/Proposals:Orientation#DICOM_LPS_Orientation_-_Full_Body_Graphical_Example

and select the position of the slice origin accordingly.

For example, for a Coronal Slice, you want the
slice origin to be at:

X := Ox
Y := (will determine the slice position)
Z := Oz + Sz * Dz


Where (Ox, Oy, Oz ) are the original coordinates of
an Axial Scan (e.g. the LPS coordinate system).

(Sx,Sy,Sz) is the number of pixels along X,Y,Z

and

(Dx,Dy,Dz) is the pixel spacing along X,Y,Z


Finally, if you are only interested in the orientation
change, and not intend to resample the image
(change its resolution), then you should make sure
that you also permute the pixel spacing according
the original values of (Dx,Dy,Dz).


Regards,


     Luis



-----------------------------------------------------------
On Mon, May 17, 2010 at 11:50 AM, Daanen Vincent <daanen at koelis.com> wrote:
> Dear Itk users,
>
> I'm facing to a problem to use ItkImageResampler to extract slice from a  3D
> image.
> The goal is to extract "AXIAL", "CORONAL" or "SAGITTAL" slice from a 3D
> image, whatever is its native orientation.
>
> I succeed in extracting Axial, coronal or sagittal slice from a Axial 3D
> image but I can't with a coronal 3D image.
>
> Here's what I'm doing.
>
> Say I have a 3D image which is coronal (or closely). Given my input image
> orientation close to
> X_Img_3D = (1,0,0), Y_Img_3D = (0,0,-1), Z_Img_3D=(0,1,0), when I want to
> extract a patient "axial" image, I have to set output image orientation
> to
>  * X_Slice = (1,0,0) ( -> = X_Img_3D)
>  * Y_Slice = (0,1,0) ( -> = Z_Img_3D)
>  * Z_Slice = (0,0,1) ( -> = -Y_Img_3D)
>
>
> // Get first "Axial"  image
> start[0] = 0;           // first index on X
> start[1] = 0;           // first index on Y
> start[2] = 0;           // first index on Z
>
> itk::Point<double,3> l_Org;
> m_pImg->TransformIndexToPhysicalPoint(start, l_Org);            // m_pImg is
> the input 3D image
> m_pResampler->SetOutputOrigin(l_Org);
>
> //! output image Direction
> TItkInputImage::DirectionType l_InputImgDir = m_pImg->GetDirection();
> TItkInputImage::DirectionType l_SliceDir;
>
> // vector X of slice is vector X of 3D image
> l_SliceDir(0,0) = l_Inputdir(0,0);
> l_SliceDir(1,0) = l_Inputdir(1,0);
> l_SliceDir(2,0) = l_Inputdir(2,0);
>
> // vector Y of slice is vector Z of 3D image
> l_SliceDir(0,1) = l_Inputdir(0,2);
> l_SliceDir(1,1) = l_Inputdir(1,2);
> l_SliceDir(2,1) = l_Inputdir(2,2);
>
> // vector Z of slice is vector -Y of 3D image
> l_SliceDir(0,2) = -l_Inputdir(0,1);
> l_SliceDir(1,2) = -l_Inputdir(1,1);
> l_SliceDir(2,2) = -l_Inputdir(2,1);
> m_pResampler->SetOutputDirection(l_SliceDir);
>
> // Spacing -> use minimun scale to get a isotropic slice.
> const TItkInputImage::SpacingType &l_InSpacing = m_pImg->GetSpacing();
> double l_IsoSpacing =
> __min(l_InSpacing[0],__min(l_InSpacing[1],l_InSpacing[2]));
>
> TItkOutputImage::SpacingType l_OutSpacing;
> l_OutSpacing[0] = l_IsoSpacing;
> l_OutSpacing[1] = l_IsoSpacing;
> l_OutSpacing[2] = l_IsoSpacing;
> m_pResampler->SetOutputSpacing(l_OutSpacing);
>
> // output image size
> TItkOutputImage::SizeType l_OutputImgSize = inputRegion.GetSize();
> TItkInputImage::SizeType size = inputRegion.GetSize();
>
> l_OutputImgSize[0] = size[0];                                           //
> width of axial image  = width of native image
> l_OutputImgSize[1] = size[2]*l_InSpacing[1]/l_IsoSpacing;               //
> height of axial image = Depth of native image * slicethickness/pixelsize
> l_OutputImgSize[2] = 1; // size along Z
>
> m_pResampler->SetSize( l_OutputImgSize );
>
>
>
> Doing so, I get an image with a small part of the native image :(.
> I must be missing something but i can't see what;
>
> Is there a itk::ImageResampleFileter guru who could help me ?
>
> Thanks
>
>  Vince
>
>
> --------------------------------------------
> Vincent Daanen, PhD
> D&D Manager
>
> --------------------------------------------
> KOELIS
> 5, avenue du Grand Sablon 38700 La Tronche
> www.koelis.com  -  daanen at koelis.com
> Tel .+33(0) 476637588 Fax .+33(0) 476637592
> --------------------------------------------
>
> CONFIDENTIALITY This e-mail and any attachments are confidential and may
> also be privileged. If you are not the named recipient, please notify the
> sender immediately and do not disclose the contents to another person, use
> it for any purpose, or store or copy the information in any medium.
>
> "Les problèmes ne peuvent être résolus par ceux dont l'horizon se limite aux
> réalités quotidiennes,  mais par ceux qui rêvent de choses qui n'ont jamais
> existé et qui se disent : Pourquoi Pas ?" (J-F Kennedy, 1963).
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>