Yeah, it's probably the direction cosines. One thing that one has to keep in mind when using the B-spline filter is that the output is an image in terms of implementation output. However, that image, for algorithmic purposes, is defined in terms of its parametric coordinates which really have no notion of directionality in physical space. So although one can set the directions for the output of the filter, internally the filter does nothing with that information. It simply passes it to the output. It's not a problem if you keep this in mind and handle it appropriately like, for example, the way you're doing by resetting the directions and then performing your B-spline computations. Let me know how it goes. <div>
<br><br><div class="gmail_quote">On Wed, Apr 7, 2010 at 4:01 PM, Andriy Fedorov <span dir="ltr"><<a href="mailto:fedorov@bwh.harvard.edu">fedorov@bwh.harvard.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
Hi Nick,<br>
<br>
I am debugging this right now.<br>
<br>
All of the points that I add to the filter are within the output<br>
domain, I do have an explicit check for this.<br>
<br>
Regarding the two suggestions you made:<br>
<br>
1) my data indeed has directional cosines that are not identity. Does<br>
this pose a problem for the filter? In principle, it is not essential<br>
to keep non-identity for the output volume, since all voxels will be<br>
interpolated anyway. Maybe I should try resetting directions for the<br>
output, if everything fails.<br>
<br>
2) I actually checked the index that corresponds to the point that<br>
causes the error, it is [265, 197, 241], while the output size is<br>
[286, 241, 274], so it's not on the boundary. Perhaps directions is<br>
indeed the root of the problem...<br>
<br>
I will update the list if I have any significant progress.<br>
<br>
Thanks<br>
<br>
--<br>
Andriy Fedorov, Ph.D.<br>
<br>
Research Fellow<br>
Brigham and Women's Hospital<br>
Harvard Medical School<br>
75 Francis Street<br>
Boston, MA 02115 USA<br>
<a href="mailto:fedorov@bwh.harvard.edu">fedorov@bwh.harvard.edu</a><br>
<div><div></div><div class="h5"><br>
<br>
<br>
On Wed, Apr 7, 2010 at 15:50, Nick Tustison <<a href="mailto:ntustison@wustl.edu">ntustison@wustl.edu</a>> wrote:<br>
> It looks like the size of your output image is too small. When you set up<br>
> the image domain (size, spacing, etc.) in the B-spline filter, you need to<br>
> make sure that all of the points that constitute the input to that filter<br>
> are within that domain. Some of the related bugs that have bitten me in the<br>
> past are 1) dealing with image data with non-identity directional cosines<br>
> and 2) round-off errors such that the point is right on one of the<br>
> boundaries. Let me know if that helps.<br>
><br>
> On Wed, Apr 7, 2010 at 2:03 PM, Andriy Fedorov <<a href="mailto:fedorov@bwh.harvard.edu">fedorov@bwh.harvard.edu</a>><br>
> wrote:<br>
>><br>
>> Nick,<br>
>><br>
>> I am getting the following error trying to implement this approach:<br>
>><br>
>> itk::ERROR: PointSetToImageFilter(0x39c0dc0): The reparameterized<br>
>> point component 1.00002 is outside the corresponding parametric domain<br>
>> of [0, 1].<br>
>><br>
>> Here's the relevant part of my code:<br>
>><br>
>> ----><br>
>><br>
>> ConstIterType it1(image1, image1->GetBufferedRegion());<br>
>><br>
>> PointSetType::Pointer pointSet = PointSetType::New();<br>
>> unsigned long numPoints = 0;<br>
>><br>
>> InputImageType::SizeType imageSize =<br>
>> image1->GetBufferedRegion().GetSize();<br>
>> std::cout << "Processing image1. Size: " << imageSize << std::endl;<br>
>><br>
>> for(it1.GoToBegin();!it1.IsAtEnd();++it1){<br>
>> PointSetType::PointType pt;<br>
>> InputImageType::PointType ipt;<br>
>> InputImageType::IndexType idx;<br>
>> PointDataType ptData;<br>
>><br>
>> image1->TransformIndexToPhysicalPoint(it1.GetIndex(), ipt);<br>
>> if(!output->TransformPhysicalPointToIndex(ipt, idx))<br>
>> continue;<br>
>><br>
>> pt[0] = ipt[0];<br>
>> pt[1] = ipt[1];<br>
>> pt[2] = ipt[2];<br>
>><br>
>> ptData[0] = it1.Get();<br>
>><br>
>> pointSet->SetPoint(numPoints, pt);<br>
>> pointSet->SetPointData(numPoints, ptData);<br>
>> numPoints++;<br>
>> }<br>
>><br>
>> std::cout << "Total number of points: " << numPoints << std::endl;<br>
>><br>
>> std::cout << "Output image: " << output << std::endl;<br>
>><br>
>> SDAFilterType::Pointer sda = SDAFilterType::New();<br>
>> SDAFilterType::ArrayType ncps;<br>
>> ncps.Fill(4);<br>
>> sda->SetSplineOrder(3);<br>
>> sda->SetNumberOfControlPoints(ncps);<br>
>> sda->SetNumberOfLevels(1);<br>
>> sda->SetOrigin(output->GetOrigin());<br>
>> sda->SetSpacing(output->GetSpacing());<br>
>> sda->SetDirection(output->GetDirection());<br>
>> sda->SetSize(output->GetBufferedRegion().GetSize());<br>
>> sda->SetInput(pointSet);<br>
>> sda->Update();<br>
>><br>
>> <----<br>
>><br>
>> Any suggestions what I may be doing wrong to cause this error? I would<br>
>> appreciate any hints.<br>
>><br>
>> Andriy Fedorov<br>
>><br>
>><br>
>> On Tue, Apr 6, 2010 at 17:59, Nicholas Tustison <<a href="mailto:ntustison@gmail.com">ntustison@gmail.com</a>><br>
>> wrote:<br>
>> > You're right about the implementation although, if you had segmentation<br>
>> > masks for each input image you could trim the input by only using the voxels<br>
>> > within the mask. Since you are implementing it anyway, make sure that your<br>
>> > point set ordering consists of all the points from one image followed by all<br>
>> > the points from the second image, etc. The processing will go faster that<br>
>> > way. Let me know if you have any additional questions.<br>
>> ><br>
>> ><br>
>> > On Apr 6, 2010, at 4:48 PM, Andriy Fedorov wrote:<br>
>> ><br>
>> >> On Tue, Apr 6, 2010 at 14:45, Nick Tustison <<a href="mailto:ntustison@wustl.edu">ntustison@wustl.edu</a>><br>
>> >> wrote:<br>
>> >>> Hi Andriy,<br>
>> >>> Have you thought about using the<br>
>> >>> itkBSplineScatteredDataPointSetToImageFilter?<br>
>> >><br>
>> >> Hi Nick,<br>
>> >><br>
>> >> Yes, I did think about this.<br>
>> >><br>
>> >> However, the volumes are quite large, and, if I understand correctly<br>
>> >> how this will be implemented, I would need to add voxel<br>
>> >> coordinates/value pair for each voxel in each of the image to the<br>
>> >> point set to initialize the bspline filter. Please correct me if I am<br>
>> >> wrong. I was wondering, if there is an approach, which would just take<br>
>> >> an arbitrary number of images, and would do interpolation by reading<br>
>> >> directly from those images.<br>
>> >><br>
>> >> Perhaps such filter can be derived from the functionality already<br>
>> >> present in itkBSplineScatteredDataPointSetToImageFilter.<br>
>> >><br>
>> >> In any case, I am working on implementing the approach based on your<br>
>> >> filter right now. Let's see how this goes.<br>
>> >><br>
>> >> Thanks for the reply!<br>
>> >><br>
>> >> AF<br>
>> >><br>
>> >><br>
>> >><br>
>> >>> Nick<br>
>> >>> On Tue, Apr 6, 2010 at 1:40 PM, Andriy Fedorov<br>
>> >>> <<a href="mailto:fedorov@bwh.harvard.edu">fedorov@bwh.harvard.edu</a>><br>
>> >>> wrote:<br>
>> >>>><br>
>> >>>> Hi,<br>
>> >>>><br>
>> >>>> I have several MRI volumes of the same anatomy, which have been<br>
>> >>>> spatially aligned. These volumes have been acquired using different<br>
>> >>>> (orthogonal) acquisition direction, so they have thick slices, but in<br>
>> >>>> orthogonal planes with respect to each other. I would like to<br>
>> >>>> construct one volume, which would essentially combine the input<br>
>> >>>> volumes, interpolate over the multiple input volumes to produce an<br>
>> >>>> image of better resolution.<br>
>> >>>><br>
>> >>>> Is there something in ITK that could be used to help me in solving<br>
>> >>>> this problem? It seems I need to do some sort of weighted averaging<br>
>> >>>> over the neighboring voxels from the input volumes. Can anyone<br>
>> >>>> suggest<br>
>> >>>> how to approach this problem?<br>
>> >>>><br>
>> >>>> Thanks<br>
>> >>>><br>
>> >>>> Andriy Fedorov<br>
>> >>>> _____________________________________<br>
>> >>>> Powered by <a href="http://www.kitware.com" 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" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
>> >>>><br>
>> >>>> Kitware offers ITK Training Courses, for more information visit:<br>
>> >>>> <a href="http://www.kitware.com/products/protraining.html" target="_blank">http://www.kitware.com/products/protraining.html</a><br>
>> >>>><br>
>> >>>> Please keep messages on-topic and check the ITK FAQ at:<br>
>> >>>> <a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
>> >>>><br>
>> >>>> Follow this link to subscribe/unsubscribe:<br>
>> >>>> <a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
>> >>><br>
>> >>><br>
>> ><br>
>> ><br>
><br>
><br>
</div></div></blockquote></div><br></div>