[ITK-users] Rotation of anisotropic 3D image

[g2]

Hi all,

I am trying to rotate in 3D with a arbitrary rotation matrix a 3D image with
non uniform spacing (0.15, 0.15, 0.19). I am currently using
itk::AffineTransform because at some point in the future I would like to
introduce a translation. My question is about the spacing of the output
image. Since my rotation matrix can be anything, I cannot just re-use the
spacing of my original image. I actually do not care to much about the
output spacing, it can be the same, it can be isotropic, whatever as long as
the data makes sense.

here is a bit of code of what I've done so far:

typedef itk::Image<short, 3> Image3d;
Image3d::Pointer itkImage = getImageSomehow();
typedef itk::AffineTransform<double, 3> TransformType;
TransformType::Pointer Rt = TransformType::New();
TransformType::ParametersType params(12);
for (int i = 0; i < 9; i++){
params[i] = R[i/3][i%3]; // R is the rotation matrix of type
double[3][3]
}
params[9] = 0;
params[10] = 0;
params[11] = 0;

Rt->SetParameters(params);

typedef itk::ResampleImageFilter<Image3d, Image3d > FilterType;
FilterType::Pointer filter = FilterType::New();
typedef itk::NearestNeighborInterpolateImageFunction<Image3d, double >
InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
filter->SetInterpolator(interpolator);
filter->SetDefaultPixelValue(255);
filter->SetOutputOrigin(itkImage->GetOrigin());
filter->SetOutputSpacing(itkImage->GetSpacing());
//double outSpacing[3] = { 0.1, 0.1, 0.1 };
//filter->SetOutputSpacing(outSpacing);
filter->SetSize(itkImage->GetLargestPossibleRegion().GetSize());
filter->SetOutputDirection(itkImage->GetDirection());
filter->SetInput(itkImage);
filter->SetTransform(Rt);
filter->Update();


After this when I save filter->GetOutput() and open it with IKT-Snap I can
see my new image properly rotated. But it has the same spacing as the input
image, as specified with filter->SetOutputSpacing(itkImage->GetSpacing());
and to me this doesn't make sense. The axes are rotated and so should be the
spacing. When I try to put some other values, such as (0.1, 0.1, 0.1). The
output image is all black. I'm confused because to me those values are not
more erroneous than a plain copy of the input spacing.

Questions :
I want to rotate a generic 3D image with anisotropic spacing with a generic
rotation matrix (i.e. not around a single axis and not with a n*90° angle)
so that any physical point P in the original volumes maps to R*P in the
final one. How should I proceed ?

thanks



--
Sent from: http://itk-insight-users.2283740.n2.nabble.com/
_____________________________________
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.php

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

[Francois Budin]

Hello Antoine,

What you are doing is correct. What I would do is that I would use the
smallest value of the spacing and create an isotropic output image (similar
to what you are doing with setting your output spacing to {0.1,0.1,0.1},
but in an automatic manner.
The only reason I can think of why you are getting an enterely black output
image is maybe because your output image space does not cover your
transformed image. If you decrease the spacing, you also need to adjust the
image size to make sure that your transformed image will be contained in
the output image space.
I created a project a long time ago that computes the output image space
based on an input image and a transform [1]. Maybe looking at that project
will help you.

Hope this helps,
Francois

[1]
https://github.com/fbudin69500/ITKTransformTools/blob/master/GetNewSizeAndOrigin.cxx

Post by g2
Hi all,
I am trying to rotate in 3D with a arbitrary rotation matrix a 3D image with
non uniform spacing (0.15, 0.15, 0.19). I am currently using
itk::AffineTransform because at some point in the future I would like to
introduce a translation. My question is about the spacing of the output
image. Since my rotation matrix can be anything, I cannot just re-use the
spacing of my original image. I actually do not care to much about the
output spacing, it can be the same, it can be isotropic, whatever as long as
the data makes sense.
typedef itk::Image<short, 3> Image3d;
Image3d::Pointer itkImage = getImageSomehow();
typedef itk::AffineTransform<double, 3> TransformType;
TransformType::Pointer Rt = TransformType::New();
TransformType::ParametersType params(12);
for (int i = 0; i < 9; i++){
params[i] = R[i/3][i%3]; // R is the rotation matrix of type
double[3][3]
}
params[9] = 0;
params[10] = 0;
params[11] = 0;
Rt->SetParameters(params);
typedef itk::ResampleImageFilter<Image3d, Image3d > FilterType;
FilterType::Pointer filter = FilterType::New();
typedef itk::NearestNeighborInterpolateImageFunction<Image3d, double >
InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
filter->SetInterpolator(interpolator);
filter->SetDefaultPixelValue(255);
filter->SetOutputOrigin(itkImage->GetOrigin());
filter->SetOutputSpacing(itkImage->GetSpacing());
//double outSpacing[3] = { 0.1, 0.1, 0.1 };
//filter->SetOutputSpacing(outSpacing);
filter->SetSize(itkImage->GetLargestPossibleRegion().GetSize());
filter->SetOutputDirection(itkImage->GetDirection());
filter->SetInput(itkImage);
filter->SetTransform(Rt);
filter->Update();
After this when I save filter->GetOutput() and open it with IKT-Snap I can
see my new image properly rotated. But it has the same spacing as the input
image, as specified with filter->SetOutputSpacing(itkImage->GetSpacing());
and to me this doesn't make sense. The axes are rotated and so should be the
spacing. When I try to put some other values, such as (0.1, 0.1, 0.1). The
output image is all black. I'm confused because to me those values are not
more erroneous than a plain copy of the input spacing.
I want to rotate a generic 3D image with anisotropic spacing with a generic
rotation matrix (i.e. not around a single axis and not with a n*90° angle)
so that any physical point P in the original volumes maps to R*P in the
final one. How should I proceed ?
thanks
--
Sent from: http://itk-insight-users.2283740.n2.nabble.com/
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

[Antoine]

Hi Francois,

This was an issue indeed. Varying the isotropic spacing to from
{0.1,0.1,0.1} to {0.5, 0.5, 0.5} I can see the corresponding output images
"drift" as the spacing increase. That covers the black image.


My next issue is to get pixel coordinates of a 3d point P'_pix in physical
space after transformation, when all I know about that point is its
physical position P_phys before rotation.

Mathematically it's quite straightforward in physical space:
P'_phys = R*P_phys

But then I need to convert this to pixel coordinates. I tried
using filter->GetOutput()->TransformPhysicalPointToIndex(Pp, Pp_pix);
but the result unsatisfactory, the output coordinates is far from what I
read by pointing at the same point in ITK-Snap.
The result I get are not off by much but they are not close either. ie
(21,891,1027) against (319, 22, 210).

compare to the code shown above I've tried setting filter output origin and
direction not to the ones of the input image but to their rotated variant.
Bellow is some pseudo-code that illustrate this.
filter->SetOutputOrigin(R*itkImage->GetOrigin());
filter->SetOutputDirection(R*itkImage->GetDirection());


Ultimately my goal is to crop a region of 1cm around a arbitrary 3D
segment. Knowing only the physical position of the segment end points. See
my awesome drawing for a 2D example of what I try to achieve.
Loading Image...
Before cropping the image I need to re-orient it so that the segment aligns
with one of the image axis, I chose X (1,0,0) but it doesn't matter. I got
that part working, as I can see my segment aligned with the chosen axis in
the temporary output image. Also I know computing such a rotation is an
under-constrained problem and there are an infinity of solutions, but I'm
fine with any of those.

But now I struggle getting proper "start" and "size" values for this piece
of code coming afterward:

Image3d::IndexType start;
Image3d::SizeType size;

//---
// missing code to get start and size right, in pixel space.
//---

Image3d::RegionType desiredRegion(start, size);
typedef itk::ExtractImageFilter< Image3d, Image3d > CropFilterType;
CropFilterType::Pointer cropFilter = CropFilterType::New();
cropFilter->SetExtractionRegion(desiredRegion);
cropFilter->SetInput(filter->GetOutput());
cropFilter->SetDirectionCollapseToIdentity();
cropFilter->Update();
Image3d::Pointer cropOut = cropFilter->GetOutput();


when I set by hand pixel values I read from ITK-snap for start and size,
this works like a charm. But I struggle to get those values automatiocally.


Cheers,

A.

https://i.imgur.com/39l3raP.png

Post by Francois Budin
Hello Antoine,
What you are doing is correct. What I would do is that I would use the
smallest value of the spacing and create an isotropic output image (similar
to what you are doing with setting your output spacing to {0.1,0.1,0.1},
but in an automatic manner.
The only reason I can think of why you are getting an enterely black
output image is maybe because your output image space does not cover your
transformed image. If you decrease the spacing, you also need to adjust the
image size to make sure that your transformed image will be contained in
the output image space.
I created a project a long time ago that computes the output image space
based on an input image and a transform [1]. Maybe looking at that project
will help you.
Hope this helps,
Francois
[1] https://github.com/fbudin69500/ITKTransformTools/blob/master/
GetNewSizeAndOrigin.cxx

Post by g2
Hi all,
I am trying to rotate in 3D with a arbitrary rotation matrix a 3D image with
non uniform spacing (0.15, 0.15, 0.19). I am currently using
itk::AffineTransform because at some point in the future I would like to
introduce a translation. My question is about the spacing of the output
image. Since my rotation matrix can be anything, I cannot just re-use the
spacing of my original image. I actually do not care to much about the
output spacing, it can be the same, it can be isotropic, whatever as long as
the data makes sense.
typedef itk::Image<short, 3> Image3d;
Image3d::Pointer itkImage = getImageSomehow();
typedef itk::AffineTransform<double, 3> TransformType;
TransformType::Pointer Rt = TransformType::New();
TransformType::ParametersType params(12);
for (int i = 0; i < 9; i++){
params[i] = R[i/3][i%3]; // R is the rotation matrix of type
double[3][3]
}
params[9] = 0;
params[10] = 0;
params[11] = 0;
Rt->SetParameters(params);
typedef itk::ResampleImageFilter<Image3d, Image3d > FilterType;
FilterType::Pointer filter = FilterType::New();
typedef itk::NearestNeighborInterpolateImageFunction<Image3d, double >
InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
filter->SetInterpolator(interpolator);
filter->SetDefaultPixelValue(255);
filter->SetOutputOrigin(itkImage->GetOrigin());
filter->SetOutputSpacing(itkImage->GetSpacing());
//double outSpacing[3] = { 0.1, 0.1, 0.1 };
//filter->SetOutputSpacing(outSpacing);
filter->SetSize(itkImage->GetLargestPossibleRegion().GetSize());
filter->SetOutputDirection(itkImage->GetDirection());
filter->SetInput(itkImage);
filter->SetTransform(Rt);
filter->Update();
After this when I save filter->GetOutput() and open it with IKT-Snap I can
see my new image properly rotated. But it has the same spacing as the input
image, as specified with filter->SetOutputSpacing(itkIm
age->GetSpacing());
and to me this doesn't make sense. The axes are rotated and so should be the
spacing. When I try to put some other values, such as (0.1, 0.1, 0.1). The
output image is all black. I'm confused because to me those values are not
more erroneous than a plain copy of the input spacing.
I want to rotate a generic 3D image with anisotropic spacing with a generic
rotation matrix (i.e. not around a single axis and not with a n*90° angle)
so that any physical point P in the original volumes maps to R*P in the
final one. How should I proceed ?
thanks
--
Sent from: http://itk-insight-users.2283740.n2.nabble.com/
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

[Dženan Zukić]

Hi Antoine,

filter->GetOutput()->TransformPhysicalPointToIndex(Pp, Pp_pix); should
work, but it can only be called after filter->Update();

If this does not resolve your issue, can you provide a complete example
with some data?

Regards,
DÅŸenan Zukić, PhD, Senior R&D Engineer, Kitware (Carrboro, N.C.)

Post by Antoine
Hi Francois,
This was an issue indeed. Varying the isotropic spacing to from
{0.1,0.1,0.1} to {0.5, 0.5, 0.5} I can see the corresponding output images
"drift" as the spacing increase. That covers the black image.
My next issue is to get pixel coordinates of a 3d point P'_pix in physical
space after transformation, when all I know about that point is its
physical position P_phys before rotation.
P'_phys = R*P_phys
But then I need to convert this to pixel coordinates. I tried
using filter->GetOutput()->TransformPhysicalPointToIndex(Pp, Pp_pix);
but the result unsatisfactory, the output coordinates is far from what I
read by pointing at the same point in ITK-Snap.
The result I get are not off by much but they are not close either. ie
(21,891,1027) against (319, 22, 210).
compare to the code shown above I've tried setting filter output origin
and direction not to the ones of the input image but to their rotated
variant. Bellow is some pseudo-code that illustrate this.
filter->SetOutputOrigin(R*itkImage->GetOrigin());
filter->SetOutputDirection(R*itkImage->GetDirection());
Ultimately my goal is to crop a region of 1cm around a arbitrary 3D
segment. Knowing only the physical position of the segment end points. See
my awesome drawing for a 2D example of what I try to achieve.
https://i.imgur.com/39l3raP.png
Before cropping the image I need to re-orient it so that the segment
aligns with one of the image axis, I chose X (1,0,0) but it doesn't
matter. I got that part working, as I can see my segment aligned with the
chosen axis in the temporary output image. Also I know computing such a
rotation is an under-constrained problem and there are an infinity of
solutions, but I'm fine with any of those.
But now I struggle getting proper "start" and "size" values for this piece
Image3d::IndexType start;
Image3d::SizeType size;
//---
// missing code to get start and size right, in pixel space.
//---
Image3d::RegionType desiredRegion(start, size);
typedef itk::ExtractImageFilter< Image3d, Image3d > CropFilterType;
CropFilterType::Pointer cropFilter = CropFilterType::New();
cropFilter->SetExtractionRegion(desiredRegion);
cropFilter->SetInput(filter->GetOutput());
cropFilter->SetDirectionCollapseToIdentity();
cropFilter->Update();
Image3d::Pointer cropOut = cropFilter->GetOutput();
when I set by hand pixel values I read from ITK-snap for start and size,
this works like a charm. But I struggle to get those values automatiocally.
Cheers,
A.
https://i.imgur.com/39l3raP.png

Post by Francois Budin
Hello Antoine,
What you are doing is correct. What I would do is that I would use the
smallest value of the spacing and create an isotropic output image (similar
to what you are doing with setting your output spacing to {0.1,0.1,0.1},
but in an automatic manner.
The only reason I can think of why you are getting an enterely black
output image is maybe because your output image space does not cover your
transformed image. If you decrease the spacing, you also need to adjust the
image size to make sure that your transformed image will be contained in
the output image space.
I created a project a long time ago that computes the output image space
based on an input image and a transform [1]. Maybe looking at that project
will help you.
Hope this helps,
Francois
[1] https://github.com/fbudin69500/ITKTransformTools/blob/
master/GetNewSizeAndOrigin.cxx

Post by g2
Hi all,
I am trying to rotate in 3D with a arbitrary rotation matrix a 3D image with
non uniform spacing (0.15, 0.15, 0.19). I am currently using
itk::AffineTransform because at some point in the future I would like to
introduce a translation. My question is about the spacing of the output
image. Since my rotation matrix can be anything, I cannot just re-use the
spacing of my original image. I actually do not care to much about the
output spacing, it can be the same, it can be isotropic, whatever as long as
the data makes sense.
typedef itk::Image<short, 3> Image3d;
Image3d::Pointer itkImage = getImageSomehow();
typedef itk::AffineTransform<double, 3> TransformType;
TransformType::Pointer Rt = TransformType::New();
TransformType::ParametersType params(12);
for (int i = 0; i < 9; i++){
params[i] = R[i/3][i%3]; // R is the rotation matrix of type
double[3][3]
}
params[9] = 0;
params[10] = 0;
params[11] = 0;
Rt->SetParameters(params);
typedef itk::ResampleImageFilter<Image3d, Image3d > FilterType;
FilterType::Pointer filter = FilterType::New();
typedef itk::NearestNeighborInterpolateImageFunction<Image3d, double >
InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
filter->SetInterpolator(interpolator);
filter->SetDefaultPixelValue(255);
filter->SetOutputOrigin(itkImage->GetOrigin());
filter->SetOutputSpacing(itkImage->GetSpacing());
//double outSpacing[3] = { 0.1, 0.1, 0.1 };
//filter->SetOutputSpacing(outSpacing);
filter->SetSize(itkImage->GetLargestPossibleRegion().GetSize());
filter->SetOutputDirection(itkImage->GetDirection());
filter->SetInput(itkImage);
filter->SetTransform(Rt);
filter->Update();
After this when I save filter->GetOutput() and open it with IKT-Snap I can
see my new image properly rotated. But it has the same spacing as the input
image, as specified with filter->SetOutputSpacing(itkIm
age->GetSpacing());
and to me this doesn't make sense. The axes are rotated and so should be the
spacing. When I try to put some other values, such as (0.1, 0.1, 0.1). The
output image is all black. I'm confused because to me those values are not
more erroneous than a plain copy of the input spacing.
I want to rotate a generic 3D image with anisotropic spacing with a generic
rotation matrix (i.e. not around a single axis and not with a n*90° angle)
so that any physical point P in the original volumes maps to R*P in the
final one. How should I proceed ?
thanks
--
Sent from: http://itk-insight-users.2283740.n2.nabble.com/
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

[Antoine]

Hi DÅŸenan,
I do call filter->Update() before filter->GetOutput()->
TransformPhysicalPointToIndex(Pp, Pp_pix).

I'll try to get you a meaningful volume data along with a standalone piece
of code by tomorrow. (company code + anonymous medical data = cannot share
as is). And I guess it's going to help me get things straight and clean.

thanks for your interest.

A.

Post by Dženan Zukić
Hi Antoine,
filter->GetOutput()->TransformPhysicalPointToIndex(Pp, Pp_pix); should
work, but it can only be called after filter->Update();
If this does not resolve your issue, can you provide a complete example
with some data?
Regards,
DÅŸenan Zukić, PhD, Senior R&D Engineer, Kitware (Carrboro, N.C.)

Post by Antoine
Hi Francois,
This was an issue indeed. Varying the isotropic spacing to from
{0.1,0.1,0.1} to {0.5, 0.5, 0.5} I can see the corresponding output images
"drift" as the spacing increase. That covers the black image.
My next issue is to get pixel coordinates of a 3d point P'_pix in
physical space after transformation, when all I know about that point is
its physical position P_phys before rotation.
P'_phys = R*P_phys
But then I need to convert this to pixel coordinates. I tried
using filter->GetOutput()->TransformPhysicalPointToIndex(Pp, Pp_pix);
but the result unsatisfactory, the output coordinates is far from what I
read by pointing at the same point in ITK-Snap.
The result I get are not off by much but they are not close either. ie
(21,891,1027) against (319, 22, 210).
compare to the code shown above I've tried setting filter output origin
and direction not to the ones of the input image but to their rotated
variant. Bellow is some pseudo-code that illustrate this.
filter->SetOutputOrigin(R*itkImage->GetOrigin());
filter->SetOutputDirection(R*itkImage->GetDirection());
Ultimately my goal is to crop a region of 1cm around a arbitrary 3D
segment. Knowing only the physical position of the segment end points. See
my awesome drawing for a 2D example of what I try to achieve.
https://i.imgur.com/39l3raP.png
Before cropping the image I need to re-orient it so that the segment
aligns with one of the image axis, I chose X (1,0,0) but it doesn't
matter. I got that part working, as I can see my segment aligned with the
chosen axis in the temporary output image. Also I know computing such a
rotation is an under-constrained problem and there are an infinity of
solutions, but I'm fine with any of those.
But now I struggle getting proper "start" and "size" values for this
Image3d::IndexType start;
Image3d::SizeType size;
//---
// missing code to get start and size right, in pixel space.
//---
Image3d::RegionType desiredRegion(start, size);
typedef itk::ExtractImageFilter< Image3d, Image3d > CropFilterType;
CropFilterType::Pointer cropFilter = CropFilterType::New();
cropFilter->SetExtractionRegion(desiredRegion);
cropFilter->SetInput(filter->GetOutput());
cropFilter->SetDirectionCollapseToIdentity();
cropFilter->Update();
Image3d::Pointer cropOut = cropFilter->GetOutput();
when I set by hand pixel values I read from ITK-snap for start and size,
this works like a charm. But I struggle to get those values automatiocally.
Cheers,
A.
https://i.imgur.com/39l3raP.png

Post by Francois Budin
Hello Antoine,
What you are doing is correct. What I would do is that I would use the
smallest value of the spacing and create an isotropic output image (similar
to what you are doing with setting your output spacing to {0.1,0.1,0.1},
but in an automatic manner.
The only reason I can think of why you are getting an enterely black
output image is maybe because your output image space does not cover your
transformed image. If you decrease the spacing, you also need to adjust the
image size to make sure that your transformed image will be contained in
the output image space.
I created a project a long time ago that computes the output image space
based on an input image and a transform [1]. Maybe looking at that project
will help you.
Hope this helps,
Francois
[1] https://github.com/fbudin69500/ITKTransformTools/blob/master
/GetNewSizeAndOrigin.cxx

Post by g2
Hi all,
I am trying to rotate in 3D with a arbitrary rotation matrix a 3D image with
non uniform spacing (0.15, 0.15, 0.19). I am currently using
itk::AffineTransform because at some point in the future I would like to
introduce a translation. My question is about the spacing of the output
image. Since my rotation matrix can be anything, I cannot just re-use the
spacing of my original image. I actually do not care to much about the
output spacing, it can be the same, it can be isotropic, whatever as long as
the data makes sense.
typedef itk::Image<short, 3> Image3d;
Image3d::Pointer itkImage = getImageSomehow();
typedef itk::AffineTransform<double, 3> TransformType;
TransformType::Pointer Rt = TransformType::New();
TransformType::ParametersType params(12);
for (int i = 0; i < 9; i++){
params[i] = R[i/3][i%3]; // R is the rotation matrix of type
double[3][3]
}
params[9] = 0;
params[10] = 0;
params[11] = 0;
Rt->SetParameters(params);
typedef itk::ResampleImageFilter<Image3d, Image3d > FilterType;
FilterType::Pointer filter = FilterType::New();
typedef itk::NearestNeighborInterpolateImageFunction<Image3d, double >
InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
filter->SetInterpolator(interpolator);
filter->SetDefaultPixelValue(255);
filter->SetOutputOrigin(itkImage->GetOrigin());
filter->SetOutputSpacing(itkImage->GetSpacing());
//double outSpacing[3] = { 0.1, 0.1, 0.1 };
//filter->SetOutputSpacing(outSpacing);
filter->SetSize(itkImage->GetLargestPossibleRegion().GetSize());
filter->SetOutputDirection(itkImage->GetDirection());
filter->SetInput(itkImage);
filter->SetTransform(Rt);
filter->Update();
After this when I save filter->GetOutput() and open it with IKT-Snap I can
see my new image properly rotated. But it has the same spacing as the input
image, as specified with filter->SetOutputSpacing(itkIm
age->GetSpacing());
and to me this doesn't make sense. The axes are rotated and so should be the
spacing. When I try to put some other values, such as (0.1, 0.1, 0.1). The
output image is all black. I'm confused because to me those values are not
more erroneous than a plain copy of the input spacing.
I want to rotate a generic 3D image with anisotropic spacing with a generic
rotation matrix (i.e. not around a single axis and not with a n*90° angle)
so that any physical point P in the original volumes maps to R*P in the
final one. How should I proceed ?
thanks
--
Sent from: http://itk-insight-users.2283740.n2.nabble.com/
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

[Francois Budin]

Hello Antoine,

I think you are missing a step in your computation:
1) You need to trnasform your point before rotation to find its position
after transformation. You can directly use the affine transform for that:
affineTransform->TransformPoint(my_point)
2) Compute the index of the transformed point:
filter->GetOutput()->TransformPhysicalPointToIndex(Pp, Pp_pix);

I think you forgot to compute 1)

Hope this helps,
Francois

Post by Antoine
Hi DÅŸenan,
I do call filter->Update() before filter->GetOutput()->Tr
ansformPhysicalPointToIndex(Pp, Pp_pix).
I'll try to get you a meaningful volume data along with a standalone piece
of code by tomorrow. (company code + anonymous medical data = cannot share
as is). And I guess it's going to help me get things straight and clean.
thanks for your interest.
A.

Post by Dženan Zukić
Hi Antoine,
filter->GetOutput()->TransformPhysicalPointToIndex(Pp, Pp_pix); should
work, but it can only be called after filter->Update();
If this does not resolve your issue, can you provide a complete example
with some data?
Regards,
DÅŸenan Zukić, PhD, Senior R&D Engineer, Kitware (Carrboro, N.C.)

Post by Antoine
Hi Francois,
This was an issue indeed. Varying the isotropic spacing to from
{0.1,0.1,0.1} to {0.5, 0.5, 0.5} I can see the corresponding output images
"drift" as the spacing increase. That covers the black image.
My next issue is to get pixel coordinates of a 3d point P'_pix in
physical space after transformation, when all I know about that point is
its physical position P_phys before rotation.
P'_phys = R*P_phys
But then I need to convert this to pixel coordinates. I tried
using filter->GetOutput()->TransformPhysicalPointToIndex(Pp, Pp_pix);
but the result unsatisfactory, the output coordinates is far from what
I read by pointing at the same point in ITK-Snap.
The result I get are not off by much but they are not close either. ie
(21,891,1027) against (319, 22, 210).
compare to the code shown above I've tried setting filter output origin
and direction not to the ones of the input image but to their rotated
variant. Bellow is some pseudo-code that illustrate this.
filter->SetOutputOrigin(R*itkImage->GetOrigin());
filter->SetOutputDirection(R*itkImage->GetDirection());
Ultimately my goal is to crop a region of 1cm around a arbitrary 3D
segment. Knowing only the physical position of the segment end points. See
my awesome drawing for a 2D example of what I try to achieve.
https://i.imgur.com/39l3raP.png
Before cropping the image I need to re-orient it so that the segment
aligns with one of the image axis, I chose X (1,0,0) but it doesn't
matter. I got that part working, as I can see my segment aligned with the
chosen axis in the temporary output image. Also I know computing such a
rotation is an under-constrained problem and there are an infinity of
solutions, but I'm fine with any of those.
But now I struggle getting proper "start" and "size" values for this
Image3d::IndexType start;
Image3d::SizeType size;
//---
// missing code to get start and size right, in pixel space.
//---
Image3d::RegionType desiredRegion(start, size);
typedef itk::ExtractImageFilter< Image3d, Image3d > CropFilterType;
CropFilterType::Pointer cropFilter = CropFilterType::New();
cropFilter->SetExtractionRegion(desiredRegion);
cropFilter->SetInput(filter->GetOutput());
cropFilter->SetDirectionCollapseToIdentity();
cropFilter->Update();
Image3d::Pointer cropOut = cropFilter->GetOutput();
when I set by hand pixel values I read from ITK-snap for start and size,
this works like a charm. But I struggle to get those values automatiocally.
Cheers,
A.
https://i.imgur.com/39l3raP.png

Post by Francois Budin
Hello Antoine,
What you are doing is correct. What I would do is that I would use the
smallest value of the spacing and create an isotropic output image (similar
to what you are doing with setting your output spacing to {0.1,0.1,0.1},
but in an automatic manner.
The only reason I can think of why you are getting an enterely black
output image is maybe because your output image space does not cover your
transformed image. If you decrease the spacing, you also need to adjust the
image size to make sure that your transformed image will be contained in
the output image space.
I created a project a long time ago that computes the output image
space based on an input image and a transform [1]. Maybe looking at that
project will help you.
Hope this helps,
Francois
[1] https://github.com/fbudin69500/ITKTransformTools/blob/master
/GetNewSizeAndOrigin.cxx

Post by g2
Hi all,
I am trying to rotate in 3D with a arbitrary rotation matrix a 3D image with
non uniform spacing (0.15, 0.15, 0.19). I am currently using
itk::AffineTransform because at some point in the future I would like to
introduce a translation. My question is about the spacing of the output
image. Since my rotation matrix can be anything, I cannot just re-use the
spacing of my original image. I actually do not care to much about the
output spacing, it can be the same, it can be isotropic, whatever as long as
the data makes sense.
typedef itk::Image<short, 3> Image3d;
Image3d::Pointer itkImage = getImageSomehow();
typedef itk::AffineTransform<double, 3> TransformType;
TransformType::Pointer Rt = TransformType::New();
TransformType::ParametersType params(12);
for (int i = 0; i < 9; i++){
params[i] = R[i/3][i%3]; // R is the rotation matrix of type
double[3][3]
}
params[9] = 0;
params[10] = 0;
params[11] = 0;
Rt->SetParameters(params);
typedef itk::ResampleImageFilter<Image3d, Image3d > FilterType;
FilterType::Pointer filter = FilterType::New();
typedef itk::NearestNeighborInterpolateImageFunction<Image3d, double >
InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
filter->SetInterpolator(interpolator);
filter->SetDefaultPixelValue(255);
filter->SetOutputOrigin(itkImage->GetOrigin());
filter->SetOutputSpacing(itkImage->GetSpacing());
//double outSpacing[3] = { 0.1, 0.1, 0.1 };
//filter->SetOutputSpacing(outSpacing);
filter->SetSize(itkImage->GetLargestPossibleRegion().GetSize());
filter->SetOutputDirection(itkImage->GetDirection());
filter->SetInput(itkImage);
filter->SetTransform(Rt);
filter->Update();
After this when I save filter->GetOutput() and open it with IKT-Snap I can
see my new image properly rotated. But it has the same spacing as the input
image, as specified with filter->SetOutputSpacing(itkIm
age->GetSpacing());
and to me this doesn't make sense. The axes are rotated and so should be the
spacing. When I try to put some other values, such as (0.1, 0.1, 0.1). The
output image is all black. I'm confused because to me those values are not
more erroneous than a plain copy of the input spacing.
I want to rotate a generic 3D image with anisotropic spacing with a generic
rotation matrix (i.e. not around a single axis and not with a n*90° angle)
so that any physical point P in the original volumes maps to R*P in the
final one. How should I proceed ?
thanks
--
Sent from: http://itk-insight-users.2283740.n2.nabble.com/
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

[g2]

Hello Francois, Dženan,

Post by Francois Budin
1) You need to trnasform your point before rotation to find its position
affineTransform->TransformPoint(my_point)
filter->GetOutput()->TransformPhysicalPointToIndex(Pp, Pp_pix);
I think you forgot to compute 1)

I did not use the TransformPoint function but I did transform the points
using my own matrix multiplication method before using
transformPhysicalPointToIndex. And it gives the same results.



While cleaning up the code I discovered that the input 3D point I was
getting were actually not in proper physical space, they did not take image
origin into account, only pixel position * spacing. I'll try to see the
implication of that in my code and I'll come back to you.

Cheers,
A.




--
Sent from: http://itk-users.7.n7.nabble.com/
_____________________________________
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.php

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://public.kitware.com/mailm

[Antoine]

OK I got it to work by changing my approach.
As I wasn't able to get the rotation and translation working in a single
pass, I divided my algorithm in two steps :
- a translation so that my initial point is in the center of the image.
- a rotation around the center of the image with the given rotation matrix

As a pre-processing I also reset the origin of my input image to be
(0,0,0), this values is then restored afterward, but it simplified a lot of
things for me.
This solved also the issue of guessing the size of my output image when the
feature I'm interested in is far from the center. By applying a translation
first I know that it's going to stay centred after the rotation and I do
not need to worry about my area of interest ending up out of bounds. After
these two steps I can crop from the center of the image and along the X
axis (see drawing from two posts above).
Obviously two re-sampling filters instead of one is bad practice. But in my
case it's done only once so it's not an issue.

here is a final bit of code doing exactly that :


// input data :
typedef itk::Image<short, 3> Image3d;
Image3d::PointType P1 = somePoint();
Image3d::Pointer itkImage = someImage();
double R[3][3] = someRotationMatrix();
auto inSize = itkImage->GetLargestPossibleRegion().GetSize();
Image3d::PointType origin = itkImage->GetOrigin();
origin.Fill(0);
itkImage->SetOrigin(origin);

//--- translate image
typedef itk::TranslationTransform<double, 3> TranslationTransformType;
TranslationTransformType::Pointer transform_t =
TranslationTransformType::New();
TranslationTransformType::OutputVectorType translation;

Image3d::PointType offset;
offset.Fill(0);
Image3d::IndexType centerPix;
for (int i = 0; i < 3; i++) centerPix[i] = inSize[i] / 2;
Image3d::PointType centerPhysical;
itkImage->TransformIndexToPhysicalPoint(centerPix, centerPhysical);
offset = P1 - centerPhysical;
translation[0] = offset[0];
translation[1] = offset[1];
translation[2] = offset[2];
transform_t->Translate(translation);

typedef itk::ResampleImageFilter<Image3d, Image3d> ResampleImageFilterType;
ResampleImageFilterType::Pointer resampleFilter_t =
ResampleImageFilterType::New();
resampleFilter_t->SetTransform(transform_t.GetPointer());
resampleFilter_t->SetInput(itkImage);
double outSpacing[3] = { 0.3, 0.3, 0.3 };
resampleFilter_t->SetOutputSpacing(outSpacing);
for (int i = 0; i < 3; i++) inSize[i] *= itkImage->GetSpacing()[i] /
outSpacing[i];
resampleFilter_t->SetSize(inSize);
resampleFilter_t->SetDefaultPixelValue(255);
resampleFilter_t->Update();
Image3d::Pointer translatedImg = resampleFilter_t->GetOutput();



//--- rotate translated image
typedef itk::FixedCenterOfRotationAffineTransform<double, 3> TransformType;
TransformType::Pointer Rt = TransformType::New();
TransformType::ParametersType params(12);
params.Fill(0);
for (int i = 0; i < 9; i++) params[i] = R[i/3][i%3]; // R[i%3][i/3] ==
R.tranpsose , switch % and / to alternate between R and R.t

Rt->SetParameters(params);
Image3d::PointType rotcenter;
for (int i = 0; i < 3; i++) centerPix[i] = inSize[i] / 2.;
translatedImg->TransformIndexToPhysicalPoint(centerPix, rotcenter);
Rt->SetCenter(rotcenter);

typedef itk::ResampleImageFilter<Image3d, Image3d > FilterType;
FilterType::Pointer filter = FilterType::New();
typedef itk::NearestNeighborInterpolateImageFunction<Image3d, double >
InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
filter->SetInterpolator(interpolator);
filter->SetDefaultPixelValue(255);
filter->SetOutputOrigin(origin);
filter->SetOutputSpacing(outSpacing);
filter->SetSize(inSize);
filter->SetInput(translatedImg);
filter->SetTransform(Rt);
filter->Update();

Image3d::Pointer rotatedImg = filter->GetOutput();

//--- crop
Image3d::IndexType start;
Image3d::SizeType size;

start.Fill(0);
size.Fill(10);

start = centerPix;
start[1] -= 10 / rotatedImg->GetSpacing()[1]; // 1cm bellow in Y
start[2] -= 10 / rotatedImg->GetSpacing()[2]; // 1cm bellow in Z
size[0] = axis.length() / rotatedImg->GetSpacing()[0]; // axis lenght in X
size[1] = 20 / rotatedImg->GetSpacing()[1]; // 2cm in Y
size[2] = 20 / rotatedImg->GetSpacing()[2]; // 2cm in Z

Image3d::RegionType desiredRegion(start, size);
typedef itk::ExtractImageFilter< Image3d, Image3d > CropFilterType;
CropFilterType::Pointer cropFilter = CropFilterType::New();
cropFilter->SetExtractionRegion(desiredRegion);
cropFilter->SetInput(rotatedImg);
cropFilter->SetDirectionCollapseToIdentity();
cropFilter->Update();
Image3d::Pointer cropOut = cropFilter->GetOutput();
cropOut->SetRegions(cropOut->GetLargestPossibleRegion().GetSize());


Thanks again for the help,
A.

Hello Francois, DÅŸenan,

Post by Francois Budin
1) You need to trnasform your point before rotation to find its position
affineTransform->TransformPoint(my_point)
filter->GetOutput()->TransformPhysicalPointToIndex(Pp, Pp_pix);
I think you forgot to compute 1)

I did not use the TransformPoint function but I did transform the points
using my own matrix multiplication method before using
transformPhysicalPointToIndex. And it gives the same results.
While cleaning up the code I discovered that the input 3D point I was
getting were actually not in proper physical space, they did not take image
origin into account, only pixel position * spacing. I'll try to see the
implication of that in my code and I'll come back to you.
Cheers,
A.
--
Sent from: http://itk-users.7.n7.nabble.com/
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

[Francois Budin]

Thanks for sharing your final code!

Post by Antoine
OK I got it to work by changing my approach.
As I wasn't able to get the rotation and translation working in a single
- a translation so that my initial point is in the center of the image.
- a rotation around the center of the image with the given rotation matrix
As a pre-processing I also reset the origin of my input image to be
(0,0,0), this values is then restored afterward, but it simplified a lot of
things for me.
This solved also the issue of guessing the size of my output image when
the feature I'm interested in is far from the center. By applying a
translation first I know that it's going to stay centred after the rotation
and I do not need to worry about my area of interest ending up out of
bounds. After these two steps I can crop from the center of the image and
along the X axis (see drawing from two posts above).
Obviously two re-sampling filters instead of one is bad practice. But in
my case it's done only once so it's not an issue.
typedef itk::Image<short, 3> Image3d;
Image3d::PointType P1 = somePoint();
Image3d::Pointer itkImage = someImage();
double R[3][3] = someRotationMatrix();
auto inSize = itkImage->GetLargestPossibleRegion().GetSize();
Image3d::PointType origin = itkImage->GetOrigin();
origin.Fill(0);
itkImage->SetOrigin(origin);
//--- translate image
typedef itk::TranslationTransform<double, 3> TranslationTransformType;
TranslationTransformType::Pointer transform_t =
TranslationTransformType::New();
TranslationTransformType::OutputVectorType translation;
Image3d::PointType offset;
offset.Fill(0);
Image3d::IndexType centerPix;
for (int i = 0; i < 3; i++) centerPix[i] = inSize[i] / 2;
Image3d::PointType centerPhysical;
itkImage->TransformIndexToPhysicalPoint(centerPix, centerPhysical);
offset = P1 - centerPhysical;
translation[0] = offset[0];
translation[1] = offset[1];
translation[2] = offset[2];
transform_t->Translate(translation);
typedef itk::ResampleImageFilter<Image3d, Image3d>
ResampleImageFilterType;
ResampleImageFilterType::Pointer resampleFilter_t =
ResampleImageFilterType::New();
resampleFilter_t->SetTransform(transform_t.GetPointer());
resampleFilter_t->SetInput(itkImage);
double outSpacing[3] = { 0.3, 0.3, 0.3 };
resampleFilter_t->SetOutputSpacing(outSpacing);
for (int i = 0; i < 3; i++) inSize[i] *= itkImage->GetSpacing()[i] /
outSpacing[i];
resampleFilter_t->SetSize(inSize);
resampleFilter_t->SetDefaultPixelValue(255);
resampleFilter_t->Update();
Image3d::Pointer translatedImg = resampleFilter_t->GetOutput();
//--- rotate translated image
typedef itk::FixedCenterOfRotationAffineTransform<double, 3>
TransformType;
TransformType::Pointer Rt = TransformType::New();
TransformType::ParametersType params(12);
params.Fill(0);
for (int i = 0; i < 9; i++) params[i] = R[i/3][i%3]; // R[i%3][i/3] ==
R.tranpsose , switch % and / to alternate between R and R.t
Rt->SetParameters(params);
Image3d::PointType rotcenter;
for (int i = 0; i < 3; i++) centerPix[i] = inSize[i] / 2.;
translatedImg->TransformIndexToPhysicalPoint(centerPix, rotcenter);
Rt->SetCenter(rotcenter);
typedef itk::ResampleImageFilter<Image3d, Image3d > FilterType;
FilterType::Pointer filter = FilterType::New();
typedef itk::NearestNeighborInterpolateImageFunction<Image3d, double >
InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
filter->SetInterpolator(interpolator);
filter->SetDefaultPixelValue(255);
filter->SetOutputOrigin(origin);
filter->SetOutputSpacing(outSpacing);
filter->SetSize(inSize);
filter->SetInput(translatedImg);
filter->SetTransform(Rt);
filter->Update();
Image3d::Pointer rotatedImg = filter->GetOutput();
//--- crop
Image3d::IndexType start;
Image3d::SizeType size;
start.Fill(0);
size.Fill(10);
start = centerPix;
start[1] -= 10 / rotatedImg->GetSpacing()[1]; // 1cm bellow in Y
start[2] -= 10 / rotatedImg->GetSpacing()[2]; // 1cm bellow in Z
size[0] = axis.length() / rotatedImg->GetSpacing()[0]; // axis lenght in X
size[1] = 20 / rotatedImg->GetSpacing()[1]; // 2cm in Y
size[2] = 20 / rotatedImg->GetSpacing()[2]; // 2cm in Z
Image3d::RegionType desiredRegion(start, size);
typedef itk::ExtractImageFilter< Image3d, Image3d > CropFilterType;
CropFilterType::Pointer cropFilter = CropFilterType::New();
cropFilter->SetExtractionRegion(desiredRegion);
cropFilter->SetInput(rotatedImg);
cropFilter->SetDirectionCollapseToIdentity();
cropFilter->Update();
Image3d::Pointer cropOut = cropFilter->GetOutput();
cropOut->SetRegions(cropOut->GetLargestPossibleRegion().GetSize());
Thanks again for the help,
A.

Hello Francois, DÅŸenan,

Post by Francois Budin
1) You need to trnasform your point before rotation to find its position
after transformation. You can directly use the affine transform for
affineTransform->TransformPoint(my_point)
filter->GetOutput()->TransformPhysicalPointToIndex(Pp, Pp_pix);
I think you forgot to compute 1)

I did not use the TransformPoint function but I did transform the points
using my own matrix multiplication method before using
transformPhysicalPointToIndex. And it gives the same results.
While cleaning up the code I discovered that the input 3D point I was
getting were actually not in proper physical space, they did not take image
origin into account, only pixel position * spacing. I'll try to see the
implication of that in my code and I'll come back to you.
Cheers,
A.
--
Sent from: http://itk-users.7.n7.nabble.com/
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users