<div dir="ltr"><div><div><div><div><div><div><div>Hello Antoine,<br><br></div>I think you are missing a step in your computation:<br></div>1) You need to trnasform your point before rotation to find its position after transformation. You can directly use the affine transform for that:<br></div>affineTransform->TransformPoint(my_point)<br>2) Compute the index of the transformed point:<br></div> filter->GetOutput()-><wbr>TransformPhysicalPointToIndex(<wbr>Pp, Pp_pix);<br><br></div>I think you forgot to compute 1)<br><br></div>Hope this helps,<br></div>Francois<br></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Sep 12, 2017 at 11:33 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">Hi <span style="font-family:verdana,sans-serif;font-size:12.8px">Dženan, </span><div><span style="font-family:verdana,sans-serif;font-size:12.8px">I do call filter->Update() before </span><span style="font-family:monospace,monospace;font-size:12.8px">filter->GetOutput()-></span><span style="font-family:monospace,monospace;font-size:12.8px">Tr<wbr>ansformPhysicalPointToIndex(</span><span style="font-family:monospace,monospace;font-size:12.8px">Pp<wbr>, Pp_pix).</span></div><div><span style="font-family:monospace,monospace;font-size:12.8px"><br></span></div><div><font face="monospace, monospace"><span style="font-size:12.8px">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.</span></font></div><div><span style="font-family:monospace,monospace;font-size:12.8px"><br></span></div><div><span style="font-family:verdana,sans-serif;font-size:12.8px">thanks for your interest.</span><br></div><div><span style="font-family:verdana,sans-serif;font-size:12.8px"><br></span></div><div><span style="font-family:verdana,sans-serif;font-size:12.8px">A.</span></div><div><span style="font-family:monospace,monospace;font-size:12.8px"><br></span></div><div><span style="font-family:monospace,monospace;font-size:12.8px"><br></span></div></div><div class="HOEnZb"><div class="h5"><div class="gmail_extra"><br><div class="gmail_quote">2017-09-12 17:13 GMT+02:00 Dženan Zukić <span dir="ltr"><<a href="mailto:dzenanz@gmail.com" target="_blank">dzenanz@gmail.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">Hi Antoine,</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default"><font face="monospace, monospace">filter->GetOutput()->Transform<wbr>PhysicalPointToIndex(Pp, Pp_pix);</font><font face="verdana, sans-serif"> </font><font face="verdana, sans-serif">should work, but it can only be called after </font><font face="monospace, monospace">filter->Update();</font></div><div class="gmail_default"><font face="verdana, sans-serif"><div dir="ltr"><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">If this does not resolve your issue, can you provide a complete example with some data?</div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif;font-size:small">Regards,</div><div class="gmail_default"><font face="verdana, sans-serif">Dženan Zukić, PhD, Senior R&D Engineer, Kitware (Carrboro, N.C.)</font></div></div></font></div></div><div class="m_-5137838834074090400HOEnZb"><div class="m_-5137838834074090400h5"><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Sep 12, 2017 at 10:43 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">Hi Francois,<div><br></div><div>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.</div><div><br></div><div><br></div><div>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.</div><div><br></div><div>Mathematically it's quite straightforward in physical space:</div><div>P'_phys = R*P_phys</div><div><br></div><div>But then I need to convert this to pixel coordinates. I tried using filter->GetOutput()->Tra<wbr>nsformPhysicalPointToIndex(Pp, Pp_pix);</div><div>but the result  unsatisfactory, the output coordinates is far from what I read by pointing at the same point in ITK-Snap. </div><div>The result I get are not off by much but they are not close either. ie (21,891,1027) against (319, 22, 210).</div><div><br></div><div>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. </div><div>filter->SetOutputOrigin(R*itkI<wbr>mage->GetOrigin());</div><div>filter->SetOutputDirection(R*i<wbr>tkImage->GetDirection()); </div><div><br></div><div><br></div><div>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. <a href="https://i.imgur.com/39l3raP.png" target="_blank">https://i.imgur.com/3<wbr>9l3raP.png</a></div><div>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.</div><div><br></div><div>But now I struggle getting proper "start" and "size" values for this piece of code coming afterward:</div><div><br></div><div><div>Image3d::IndexType start;</div><div>Image3d::SizeType size;</div></div><div><br></div><div>//---</div><div>// missing code to get start and size right, in pixel space.</div><div>//---</div><div><br></div><div><div>Image3d::RegionType desiredRegion(start, size);</div><div>typedef itk::ExtractImageFilter< Image3d, Image3d > CropFilterType;</div><div>CropFilterType::Pointer cropFilter = CropFilterType::New();</div><div>cropFilter->SetExtractionRegio<wbr>n(desiredRegion);</div><div>cropFilter->SetInput(filter->G<wbr>etOutput());</div><div>cropFilter->SetDirectionCollap<wbr>seToIdentity();</div><div>cropFilter->Update();</div><div>Image3d::Pointer cropOut = cropFilter->GetOutput();<br></div></div><div><br></div><div><br></div><div>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.</div><div><br></div><div><br></div><div>Cheers,</div><div><br></div><div>A.</div><div><br></div><div><a href="https://i.imgur.com/39l3raP.png" target="_blank">https://i.imgur.com/39l3raP.pn<wbr>g</a><br></div><div><br></div></div><div class="m_-5137838834074090400m_-2142207867617807457HOEnZb"><div class="m_-5137838834074090400m_-2142207867617807457h5"><div class="gmail_extra"><br><div class="gmail_quote">2017-09-12 14:55 GMT+02:00 Francois Budin <span dir="ltr"><<a href="mailto:francois.budin@kitware.com" target="_blank">francois.budin@kitware.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div><div><div>Hello Antoine,<br><br></div>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.<br></div>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.<br></div>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.<br><br></div>Hope this helps,<br></div>Francois<br><br>[1] <a href="https://github.com/fbudin69500/ITKTransformTools/blob/master/GetNewSizeAndOrigin.cxx" target="_blank">https://github.com/fbudin69500<wbr>/ITKTransformTools/blob/master<wbr>/GetNewSizeAndOrigin.cxx</a><br></div><div class="m_-5137838834074090400m_-2142207867617807457m_4557523225571091634HOEnZb"><div class="m_-5137838834074090400m_-2142207867617807457m_4557523225571091634h5"><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Sep 12, 2017 at 5:41 AM, g2 <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">Hi all,<br>
<br>
I am trying to rotate in 3D with a arbitrary rotation matrix a 3D image with<br>
non uniform spacing (0.15, 0.15, 0.19). I am currently using<br>
itk::AffineTransform because at some point in the future I would like to<br>
introduce a translation. My question is about the spacing of the output<br>
image. Since my rotation matrix can be anything, I cannot just re-use the<br>
spacing of my original image. I actually do not care to much about the<br>
output spacing, it can be the same, it can be isotropic, whatever as long as<br>
the data makes sense.<br>
<br>
here is a bit of code of what I've done so far:<br>
<br>
typedef itk::Image<short, 3> Image3d;<br>
Image3d::Pointer itkImage = getImageSomehow();<br>
typedef itk::AffineTransform<double, 3> TransformType;<br>
TransformType::Pointer Rt = TransformType::New();<br>
TransformType::ParametersType params(12);<br>
for (int i = 0; i < 9; i++){<br>
        params[i] = R[i/3][i%3]; // R is the rotation matrix of type<br>
double[3][3]<br>
}<br>
params[9] = 0;<br>
params[10] = 0;<br>
params[11] = 0;<br>
<br>
Rt->SetParameters(params);<br>
<br>
typedef itk::ResampleImageFilter<Image<wbr>3d, Image3d > FilterType;<br>
FilterType::Pointer filter = FilterType::New();<br>
typedef itk::NearestNeighborInterpolat<wbr>eImageFunction<Image3d, double ><br>
InterpolatorType;<br>
InterpolatorType::Pointer interpolator = InterpolatorType::New();<br>
filter->SetInterpolator(interp<wbr>olator);<br>
filter->SetDefaultPixelValue(2<wbr>55);<br>
filter->SetOutputOrigin(itkIma<wbr>ge->GetOrigin());<br>
filter->SetOutputSpacing(itkIm<wbr>age->GetSpacing());<br>
//double outSpacing[3] = { 0.1, 0.1, 0.1 };<br>
//filter->SetOutputSpacing(out<wbr>Spacing);<br>
filter->SetSize(itkImage->GetL<wbr>argestPossibleRegion().GetSize<wbr>());<br>
filter->SetOutputDirection(itk<wbr>Image->GetDirection());<br>
filter->SetInput(itkImage);<br>
filter->SetTransform(Rt);<br>
filter->Update();<br>
<br>
<br>
After this when I save filter->GetOutput() and open it with IKT-Snap I can<br>
see my new image properly rotated. But it has the same spacing as the input<br>
image, as specified with filter->SetOutputSpacing(itkIm<wbr>age->GetSpacing());<br>
and to me this doesn't make sense. The axes are rotated and so should be the<br>
spacing. When I try to put some other values, such as (0.1, 0.1, 0.1). The<br>
output image is all black. I'm confused because to me those values are not<br>
more erroneous than a plain copy of the input spacing.<br>
<br>
Questions :<br>
I want to rotate a generic 3D image with anisotropic spacing with a generic<br>
rotation matrix (i.e. not around a single axis and not with a n*90° angle)<br>
so that any physical point P in the original volumes maps to R*P in the<br>
final one. How should I proceed ?<br>
<br>
thanks<br>
<br>
<br>
<br>
--<br>
Sent from: <a href="http://itk-insight-users.2283740.n2.nabble.com/" rel="noreferrer" target="_blank">http://itk-insight-users.22837<wbr>40.n2.nabble.com/</a><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/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>
</blockquote></div><br></div>
</div></div></blockquote></div><br></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/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>
<br></blockquote></div><br></div>
</div></div></blockquote></div><br></div>
</div></div></blockquote></div><br></div>