[ITK-users] Rotation of anisotropic 3D image

[Antoine]

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 <francois.budin at kitware.com>:

> 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 <antoine.letouzey at gmail.com> 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
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/insight-users/attachments/20170912/e1433439/attachment.html>