Re: [ITK-users] Rotation of anisotropic 3D image

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

On Tue, Sep 12, 2017 at 10:43 AM, Antoine <[email protected]>
wrote:

> 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.
> 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
> 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
>
>
> 2017-09-12 14:55 GMT+02:00 Francois Budin <[email protected]>:
>
>> 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
>>
>> On Tue, Sep 12, 2017 at 5:41 AM, g2 <[email protected]> wrote:
>>
>>> 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(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.
>>>
>>> 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://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
>
>

_____________________________________
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