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