[Insight-users] Interpolating over multiple images

[Andriy Fedorov]

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
>>>
>>>
>
>