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

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

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

> 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.
>
>
>
> 2017-09-12 17:13 GMT+02:00 Dženan Zukić <[email protected]>:
>
>> 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