[ITK-users] Rotation of anisotropic 3D image

[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.


2017-09-13 13:16 GMT+02:00 g2 <antoine.letouzey at gmail.com>:

> Hello Francois, Dženan,
>
>
> Francois Budin-3 wrote
> > 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)
>
> 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/mailman/listinfo/insight-users
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/insight-users/attachments/20170915/8c47e07b/attachment.html>