[ITK-users] Rotation of anisotropic 3D image

[Francois Budin]

Thanks for sharing your final code!

On Fri, Sep 15, 2017 at 7:14 AM, Antoine <antoine.letouzey at gmail.com> wrote:

> 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
>>
>
>
> _____________________________________
> 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/a88e6b01/attachment.html>