[Insight-users] Interpolating over multiple images

Andriy Fedorov fedorov at bwh.harvard.edu
Fri Apr 23 10:33:20 EDT 2010 

Hi Nick and all ITK users,

I am sorry for not updating this thread for a while.

I was able to resolve the error I reported earlier by using the output
image with identity direction, and fixing another error in my code. I
have the code running now, but I have several problems. First, let me
describe the setup and approach I am using in more detail.

I have saggital, coronal and axial images (one each) of the same
anatomy, each with slice thickness about 2 mm and isotropic in-slice
resolution of about 0.2mm. I would like to create a new image with
isotropic resolution that covers the region where all three input
images overlap.

In order to set up bspline interpolation, I add each of the voxels
from the input images that lie within the output region to the set of
points that I pass to the bspline scatter filter. I set the control
point grid to be the size of the output image, bspline order to 3.

First, I tried to use isotropic spacing of 0.2 mm for the output
image. This takes a very long time to process, and in the end I am
getting something like a "grid" of interleaving planes from the input
images. I think this makes sense, because what I am trying to do is to
get values at the voxels marked "O" from the "X" voxels, as in 2d
example below, and the bspline order is 3, so it doesn't even reach
the initialized points.

XXXXXXXXXX
XOOOOOOOX
XOOOOOOOX
XOOOOOOOX
XOOOOOOOX
XOOOOOOOX
XOOOOOOOX
XOOOOOOOX
XXXXXXXXXX

I tried to reduce the resolution of the output image, but it seems
like I am getting a visually better result by simple averaging of the
intensities of the input images.

Is the approach with scatter bspline interpolation appropriate for the
problem I am trying to resolve? Maybe I am initializing the filters
incorrectly?

Are there any other solutions I could explore that are available within ITK?

I would appreciate any suggestions

AF




On Wed, Apr 7, 2010 at 16:14, Nick Tustison <ntustison at wustl.edu> wrote:
> 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
>> >> >>>
>> >> >>>
>> >> >
>> >> >
>> >
>> >
>
>

More information about the Insight-users mailing list