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

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.



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