[prev in list] [next in list] [prev in thread] [next in thread]
List: insight-users
Subject: Re: [ITK-users] Rotation of anisotropic 3D image
From: Francois Budin <francois.budin () kitware ! com>
Date: 2017-09-15 13:06:01
Message-ID: CAHwLD2Wu3RCON7kPCQdTzETrC-_gmFqUA_viU_55iTfv0tBWmw () mail ! gmail ! com
[Download RAW message or body]
[Attachment #2 (multipart/alternative)]
Thanks for sharing your final code!
On Fri, Sep 15, 2017 at 7:14 AM, Antoine <antoine.letouzey@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@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
>
>
[Attachment #5 (text/html)]
<div dir="ltr">Thanks for sharing your final code!<br></div><div \
class="gmail_extra"><br><div class="gmail_quote">On Fri, Sep 15, 2017 at 7:14 AM, \
Antoine <span dir="ltr"><<a href="mailto:antoine.letouzey@gmail.com" \
target="_blank">antoine.letouzey@gmail.com</a>></span> wrote:<br><blockquote \
class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc \
solid;padding-left:1ex"><div dir="ltr">OK I got it to work by changing my \
approach.<div>As I wasn't able to get the rotation and translation working in a \
single pass, I divided my algorithm in two steps : </div><div>- a translation so \
that my initial point is in the center of the image.</div><div>- a rotation around \
the center of the image with the given rotation matrix</div><div><br></div><div>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.</div><div>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).</div><div>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.</div><div><br></div><div>here is a final bit of code doing exactly that \
:</div><div><br></div><div><br></div><div>// input data :</div><span \
class=""><div>typedef itk::Image<short, 3> \
Image3d;<br></div></span><div><div>Image3d::PointType P1 = \
somePoint();</div></div><div>Image3d::Pointer itkImage = \
someImage();<br></div><div>double R[3][3] = someRotationMatrix();<br></div><div>auto \
inSize = itkImage-><wbr>GetLargestPossibleRegion().<wbr>GetSize();<br></div><div><div>Image3d::PointType \
origin = itkImage->GetOrigin();</div><div>origin.Fill(0);</div><div>itkImage->SetOrigin(origin);</div></div><div><div><br></div><div>//--- \
translate image</div><div><span style="white-space:pre-wrap"> </span>typedef \
itk::TranslationTransform<<wbr>double, 3> \
TranslationTransformType;</div><div><span \
style="white-space:pre-wrap"> </span>TranslationTransformType::<wbr>Pointer \
transform_t = TranslationTransformType::New(<wbr>);</div><div><span \
style="white-space:pre-wrap"> </span>TranslationTransformType::<wbr>OutputVectorType \
translation;</div><div><br></div><div><span \
style="white-space:pre-wrap"> </span>Image3d::PointType offset;</div><div><span \
style="white-space:pre-wrap"> </span>offset.Fill(0);</div><div><span \
style="white-space:pre-wrap"> </span>Image3d::IndexType centerPix;</div><div><span \
style="white-space:pre-wrap"> </span>for (int i = 0; i < 3; i++) centerPix[i] = \
inSize[i] / 2;</div><div><span \
style="white-space:pre-wrap"> </span>Image3d::PointType \
centerPhysical;</div><div><span \
style="white-space:pre-wrap"> </span>itkImage-><wbr>TransformIndexToPhysicalPoint(<wbr>centerPix, \
centerPhysical);</div><div><span style="white-space:pre-wrap"> </span>offset = P1 - \
centerPhysical;</div><div><span style="white-space:pre-wrap"> </span>translation[0] = \
offset[0];</div><div><span style="white-space:pre-wrap"> </span>translation[1] = \
offset[1];</div><div><span style="white-space:pre-wrap"> </span>translation[2] = \
offset[2];</div><div><span \
style="white-space:pre-wrap"> </span>transform_t->Translate(<wbr>translation);</div><div><br></div><div><span \
style="white-space:pre-wrap"> </span>typedef \
itk::ResampleImageFilter<<wbr>Image3d, Image3d> \
ResampleImageFilterType;</div><div><span \
style="white-space:pre-wrap"> </span>ResampleImageFilterType::<wbr>Pointer \
resampleFilter_t = ResampleImageFilterType::New()<wbr>;</div><div><span \
style="white-space:pre-wrap"> </span>resampleFilter_t-><wbr>SetTransform(transform_t.<wbr>GetPointer());</div><div><span \
style="white-space:pre-wrap"> </span>resampleFilter_t->SetInput(<wbr>itkImage);</div><div><span \
style="white-space:pre-wrap"> </span>double outSpacing[3] = { 0.3, 0.3, 0.3 \
};</div><div><span style="white-space:pre-wrap"> </span>resampleFilter_t-><wbr>SetOutputSpacing(outSpacing);</div><div><span \
style="white-space:pre-wrap"> </span>for (int i = 0; i < 3; i++) inSize[i] *= \
itkImage->GetSpacing()[i] / outSpacing[i];</div><div><span \
style="white-space:pre-wrap"> </span>resampleFilter_t->SetSize(<wbr>inSize);</div><div><span \
style="white-space:pre-wrap"> </span>resampleFilter_t-><wbr>SetDefaultPixelValue(255);</div><div><span \
style="white-space:pre-wrap"> </span>resampleFilter_t->Update();</div><div><span \
style="white-space:pre-wrap"> </span>Image3d::Pointer translatedImg = \
resampleFilter_t->GetOutput();</div><div><br></div><div><br></div><div><br></div><div><span \
style="white-space:pre-wrap"> </span>//--- rotate translated image</div><div><span \
style="white-space:pre-wrap"> </span>typedef \
itk::<wbr>FixedCenterOfRotationAffineTra<wbr>nsform<double, 3> \
TransformType;</div><span class=""><div><span \
style="white-space:pre-wrap"> </span>TransformType::Pointer Rt = \
TransformType::New();</div><div><span \
style="white-space:pre-wrap"> </span>TransformType::ParametersType \
params(12);</div></span><div><span \
style="white-space:pre-wrap"> </span>params.Fill(0);</div><div><span \
style="white-space:pre-wrap"> </span>for (int i = 0; i < 9; i++)<span \
style="white-space:pre-wrap"> </span>params[i] = R[i/3][i%3]; // R[i%3][i/3] == \
R.tranpsose , switch % and / to alternate between R and \
R.t</div><div><br></div><div><span \
style="white-space:pre-wrap"> </span>Rt->SetParameters(params);</div><div><span \
style="white-space:pre-wrap"> </span>Image3d::PointType rotcenter;</div><div><span \
style="white-space:pre-wrap"> </span>for (int i = 0; i < 3; i++) centerPix[i] = \
inSize[i] / 2.;</div><div><span \
style="white-space:pre-wrap"> </span>translatedImg-><wbr>TransformIndexToPhysicalPoint(<wbr>centerPix, \
rotcenter);</div><div><span \
style="white-space:pre-wrap"> </span>Rt->SetCenter(rotcenter);</div><span \
class=""><div><br></div><div><span style="white-space:pre-wrap"> </span>typedef \
itk::ResampleImageFilter<<wbr>Image3d, Image3d > FilterType;</div><div><span \
style="white-space:pre-wrap"> </span>FilterType::Pointer filter = \
FilterType::New();</div><div><span style="white-space:pre-wrap"> </span>typedef<span \
style="white-space:pre-wrap"> </span>itk::<wbr>NearestNeighborInterpolateImag<wbr>eFunction<Image3d, \
double > InterpolatorType;</div><div><span \
style="white-space:pre-wrap"> </span>InterpolatorType::Pointer interpolator = \
InterpolatorType::New();</div><div><span \
style="white-space:pre-wrap"> </span>filter->SetInterpolator(<wbr>interpolator);</div><div><span \
style="white-space:pre-wrap"> </span>filter->SetDefaultPixelValue(<wbr>255);</div></span><div><span \
style="white-space:pre-wrap"> </span>filter->SetOutputOrigin(<wbr>origin);</div><div><span \
style="white-space:pre-wrap"> </span>filter->SetOutputSpacing(<wbr>outSpacing);</div><div><span \
style="white-space:pre-wrap"> </span>filter->SetSize(inSize);</div><div><span \
style="white-space:pre-wrap"> </span>filter->SetInput(<wbr>translatedImg);</div><div><span \
style="white-space:pre-wrap"> </span>filter->SetTransform(Rt);</div><div><span \
style="white-space:pre-wrap"> </span>filter->Update();</div><div><br></div><div><span \
style="white-space:pre-wrap"> </span>Image3d::Pointer rotatedImg = \
filter->GetOutput();</div><div><br></div><div><span \
style="white-space:pre-wrap"> </span>//--- crop </div><span class=""><div><span \
style="white-space:pre-wrap"> </span>Image3d::IndexType start;</div><div><span \
style="white-space:pre-wrap"> </span>Image3d::SizeType \
size;</div><div><br></div></span><div><span \
style="white-space:pre-wrap"> </span>start.Fill(0);</div><div><span \
style="white-space:pre-wrap"> </span>size.Fill(10);</div><div><br></div><div><span \
style="white-space:pre-wrap"> </span>start = centerPix;</div><div><span \
style="white-space:pre-wrap"> </span>start[1] -= 10 / rotatedImg->GetSpacing()[1]; \
// 1cm bellow in Y</div><div><span style="white-space:pre-wrap"> </span>start[2] -= \
10 / rotatedImg->GetSpacing()[2]; // 1cm bellow in Z</div><div><span \
style="white-space:pre-wrap"> </span>size[0] = axis.length() / \
rotatedImg->GetSpacing()[0]; // axis lenght in X</div><div><span \
style="white-space:pre-wrap"> </span>size[1] = 20 / rotatedImg->GetSpacing()[1]; \
// 2cm in Y</div><div><span style="white-space:pre-wrap"> </span>size[2] = 20 / \
rotatedImg->GetSpacing()[2]; // 2cm in Z</div><span class=""><div><span \
style="white-space:pre-wrap"> </span></div><div><br></div><div><span \
style="white-space:pre-wrap"> </span>Image3d::RegionType desiredRegion(start, \
size);</div><div><span style="white-space:pre-wrap"> </span>typedef \
itk::ExtractImageFilter< Image3d, Image3d > CropFilterType;</div><div><span \
style="white-space:pre-wrap"> </span>CropFilterType::Pointer cropFilter = \
CropFilterType::New();</div><div><span \
style="white-space:pre-wrap"> </span>cropFilter-><wbr>SetExtractionRegion(<wbr>desiredRegion);</div></span><div><span \
style="white-space:pre-wrap"> </span>cropFilter->SetInput(<wbr>rotatedImg);</div><span \
class=""><div><span style="white-space:pre-wrap"> </span>cropFilter-><wbr>SetDirectionCollapseToIdentity<wbr>();</div><div><span \
style="white-space:pre-wrap"> </span>cropFilter->Update();</div><div><span \
style="white-space:pre-wrap"> </span>Image3d::Pointer cropOut = \
cropFilter->GetOutput();</div></span><div><span \
style="white-space:pre-wrap"> </span>cropOut->SetRegions(cropOut-><wbr>GetLarges \
tPossibleRegion().<wbr>GetSize());</div></div><div><br></div><div><br></div><div>Thanks \
again for the help,</div><div>A.</div><div><div class="h5"><div><br></div><div \
class="gmail_extra"><br><div class="gmail_quote">2017-09-13 13:16 GMT+02:00 g2 <span \
dir="ltr"><<a href="mailto:antoine.letouzey@gmail.com" \
target="_blank">antoine.letouzey@gmail.com</a>></span>:<br><blockquote \
class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid \
rgb(204,204,204);padding-left:1ex">Hello Francois, Dženan,<br> <br>
<br>
Francois Budin-3 wrote<br>
<span>> I think you are missing a step in your computation:<br>
> 1) You need to trnasform your point before rotation to find its position<br>
> after transformation. You can directly use the affine transform for that:<br>
> affineTransform->TransformPoin<wbr>t(my_point)<br>
> 2) Compute the index of the transformed point:<br>
> filter->GetOutput()->Transform<wbr>PhysicalPointToIndex(Pp, Pp_pix);<br>
><br>
> I think you forgot to compute 1)<br>
<br>
</span>I did not use the TransformPoint function but I did transform the points<br>
using my own matrix multiplication method before using<br>
transformPhysicalPointToIndex. And it gives the same results.<br>
<br>
<br>
<br>
While cleaning up the code I discovered that the input 3D point I was<br>
getting were actually not in proper physical space, they did not take image<br>
origin into account, only pixel position * spacing. I'll try to see the<br>
implication of that in my code and I'll come back to you.<br>
<br>
Cheers,<br>
A.<br>
<br>
<br>
<br>
<br>
--<br>
Sent from: <a href="http://itk-users.7.n7.nabble.com/" rel="noreferrer" \
target="_blank">http://itk-users.7.n7.nabble.c<wbr>om/</a><br> <div \
class="m_1320149395339190654gmail-m_1902352100012120978HOEnZb"><div \
class="m_1320149395339190654gmail-m_1902352100012120978h5">______________________________<wbr>_______<br>
Powered by <a href="http://www.kitware.com" rel="noreferrer" \
target="_blank">www.kitware.com</a><br> <br>
Visit other Kitware open-source projects at<br>
<a href="http://www.kitware.com/opensource/opensource.html" rel="noreferrer" \
target="_blank">http://www.kitware.com/opensou<wbr>rce/opensource.html</a><br> <br>
Kitware offers ITK Training Courses, for more information visit:<br>
<a href="http://www.kitware.com/products/protraining.php" rel="noreferrer" \
target="_blank">http://www.kitware.com/product<wbr>s/protraining.php</a><br> <br>
Please keep messages on-topic and check the ITK FAQ at:<br>
<a href="http://www.itk.org/Wiki/ITK_FAQ" rel="noreferrer" \
target="_blank">http://www.itk.org/Wiki/ITK_FA<wbr>Q</a><br> <br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://public.kitware.com/mailman/listinfo/insight-users" rel="noreferrer" \
target="_blank">http://public.kitware.com/mail<wbr>man/listinfo/insight-users</a><br> \
</div></div></blockquote></div><br></div></div></div></div> \
<br>______________________________<wbr>_______<br> Powered by <a \
href="http://www.kitware.com" rel="noreferrer" \
target="_blank">www.kitware.com</a><br> <br>
Visit other Kitware open-source projects at<br>
<a href="http://www.kitware.com/opensource/opensource.html" rel="noreferrer" \
target="_blank">http://www.kitware.com/<wbr>opensource/opensource.html</a><br> <br>
Kitware offers ITK Training Courses, for more information visit:<br>
<a href="http://www.kitware.com/products/protraining.php" rel="noreferrer" \
target="_blank">http://www.kitware.com/<wbr>products/protraining.php</a><br> <br>
Please keep messages on-topic and check the ITK FAQ at:<br>
<a href="http://www.itk.org/Wiki/ITK_FAQ" rel="noreferrer" \
target="_blank">http://www.itk.org/Wiki/ITK_<wbr>FAQ</a><br> <br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://public.kitware.com/mailman/listinfo/insight-users" rel="noreferrer" \
target="_blank">http://public.kitware.com/<wbr>mailman/listinfo/insight-users</a><br> \
<br></blockquote></div><br></div>
_____________________________________
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
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic