[Insight-users] Addition of two DeformationFields

[Miller, James V (Research)]

Monica, 
 
I forgot something in my last message.  The calls to
interpolator->Evaluate() should be within an if block that checks that the
point is within the image, for instance
 
    // interpolate a vector from the smaller image at this physical position

    if (interpolator->IsInsideBuffer(p)
       {
       v2 = interpolator->Evaluate(p);
       }
    else
       {
       v2.Fill( 0.0 );
       }
 
 

-----Original Message-----
From: Miller, James V (Research) [mailto:millerjv at crd.ge.com]
Sent: Friday, August 13, 2004 9:57 AM
To: 'Monica Lumbardi'; 'Luis Ibanez'
Cc: 'insight-users at itk.org'
Subject: RE: [Insight-users] Addition of two DeformationFields


Monica, 
 
The issue with the AddImageFilter is that it relies on the NumericTraits
class to define a
type for accumulating the results.  For instance, if you are adding to char
images, you would use an accumulate type of short to avoid overflow during
the accumulation.
 
There is not specialization for NumericTraits<Vector<T>>.  You might try
putting one in your code (you can use one of the existing specializations to
see what traits need to be defined). But I don't think this will work given
the class heirarchy for NumericTraits.
 
We could add another algorithm to specifically add two vector images.
 
However, looking back at your original question, it sounds like your two
deformation fields are at different resolutions. So the smaller one would
need to be interpolated.  You could something like
 
typedef itk::VectorLinearInterpolateImageFunction<SmallImageType>
VectorInterpolatorType;
 
VectorInterpolatorType::Pointer interpolator =
VectorInterpolatorType::New();
interpolator->SetInputImage( smallField );
 
OutputIteratorType outIt = OutputIteratorType(cumulativeField, region);
LargeInputIteratorType inIt = LargeInputIteratorType(largeField, region);
 
VectorType v1, v2;
PointType p;
while (!inIt.IsAtEnd())
{

    largeField->TransfromIndexToPhysicalPoint( inIt.GetIndex(), p );
    v1 = inIt.Get();
 
    // interpolate a vector from the smaller image at this physical position

    v2 = interpolator->Evaluate(p);
 
   // set the output vector as the sum of the two vectors
    outIt.Set( v1 + v2 );
 
    ++inIt;
   ++outIt;
}
 
I trust that adding deformation is what you really want to do.  In
particular, you are 
looking to add two vectors at the same physical position.  This is what the
above code does. If you are looking to "concatenate" or "compose" to vector
fields, where you apply one vector field to a set of positions, then apply
another vector field, then you use inner loop code that looks like 

 
{

    largeField->TransfromIndexToPhysicalPoint( inIt.GetIndex(), p );
    v1 = inIt.Get();
 
    // to compose the vector fields, apply the first vector field to the
position
    newp = p + v1;
    
    // interpolate the vector to compose with from the smaller image
    v2 = interpolator->Evaluate(newp);
 
    // apply the second vector field to position
    newp += v2;
 
    // composite vector is the difference between the original point and the
moved point
    outIt.Set( newp - p );
 
    ++inIt;
    ++outIt;
}
 
If you want, you can put in a feature request on the ITK bug tracker for a
filter that will composite two vector fields.
 

[Jim Miller] 
 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20040813/42d79edc/attachment-0001.html