[Insight-users] Interpolating over multiple images

Wed Apr 7 16:01:49 EDT 2010

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