[Insight-users] Interpolating over multiple images

[Nick Tustison]

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/388998d6/attachment-0001.htm>