[Insight-users] Interpolating over multiple images

Wed Apr 7 16:14:25 EDT 2010

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.

On Wed, Apr 7, 2010 at 4:01 PM, Andriy Fedorov <fedorov at bwh.harvard.edu>wrote:

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