[Insight-users] How to limit rotations when registering in IT K?
Eduard Schreibmann
eschreibmann at yahoo.com
Fri Oct 8 22:44:32 EDT 2004
Hello Robert,
The LBFGSBOptimizer in ITK is constrained. You can set there limits for all variables or only one, like the one representing the rotation. You do not need to write any additional code or make modifications to ITK, just use this optimizer, it is the bounded version of LBFGS. You can find some sample code on how to use this in the testing code.
Edi
Robert Maroon <robertmaroon at yahoo.com> wrote:
Thanks for the advice James,
I've been trying to implement it using the subclassing method but even though my code compiles fine it keeps calling the original code from MeanSquaresImageToImageMetric instead of the class I derived from it. I included my header file and a couple of clips from my code to see if anyone can see what I am doing wrong. I am trying to use all the functions of the parent and just replacing GetValueAndDerivative and GetValue with my own subclass.
Thanks in advance! I would really appreciate some help!
In my main project code I call on my class as follows:
I include:
#include "myReg.h"
And then to declare the metric I do:
typedef itk::myReg<FixedImageType,MovingImageType> MetricType;
The rest of the file is basically ImageRegistration5.cxx
#ifndef __myReg_h
#define __myReg_h
#include "itkImageToImageMetric.h"
#include "itkCovariantVector.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkPoint.h"
namespace itk
{
template < class TFixedImage, class TMovingImage >
class ITK_EXPORT myReg :
public MeanSquaresImageToImageMetric< TFixedImage, TMovingImage>
{
/** Get the value for single valued optimizers. */
MeasureType GetValue( const TransformParametersType & parameters ) const;
/** Get value and derivatives for multiple valued optimizers. */
void GetValueAndDerivative( const TransformParametersType & parameters,
MeasureType& Value, DerivativeType& Derivative ) const;
protected:
myReg();
virtual ~myReg() {};
};
} // end namespace itk
#endif
Here's the TXX file:
#ifndef _myReg_txx
#define _myReg_txx
#include "myReg.h"
#include "itkImageRegionConstIteratorWithIndex.h"
namespace itk
{
/*
* Constructor
*/
template <class TFixedImage, class TMovingImage>
myReg<TFixedImage,TMovingImage>
::myReg()
{
std::cout<<"TEST 1!"<<std::endl;
itkDebugMacro("Constructor");
}
/*
* Get the match Measure
*/
template <class TFixedImage, class TMovingImage>
typename myReg<TFixedImage,TMovingImage>::MeasureType
myReg<TFixedImage,TMovingImage>
::GetValue( const TransformParametersType & parameters ) const
{
itkDebugMacro("GetValue( " << parameters << " ) ");
FixedImageConstPointer fixedImage = this->GetFixedImage();
if( !fixedImage )
{
itkExceptionMacro( << "Fixed image has not been assigned" );
}
typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;
FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
typename FixedImageType::IndexType index;
MeasureType measure = NumericTraits< MeasureType >::Zero;
m_NumberOfPixelsCounted = 0;
this->SetTransformParameters( parameters );
while(!ti.IsAtEnd())
{
index = ti.GetIndex();
typename Superclass::InputPointType inputPoint;
fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
if( m_FixedImageMask && !m_FixedImageMask->IsInside( inputPoint ) )
{
++ti;
continue;
}
typename Superclass::OutputPointType transformedPoint = m_Transform->TransformPoint( inputPoint );
if( m_MovingImageMask && !m_MovingImageMask->IsInside( transformedPoint ) )
{
++ti;
continue;
}
if( m_Interpolator->IsInsideBuffer( transformedPoint ) )
{
const RealType movingValue = m_Interpolator->Evaluate( transformedPoint );
const RealType fixedValue = ti.Get();
m_NumberOfPixelsCounted++;
const RealType diff = movingValue - fixedValue;
measure += diff * diff;
}
++ti;
}
if( !m_NumberOfPixelsCounted )
{
itkExceptionMacro(<<"All the points mapped to outside of the moving image");
}
else
{
measure /= m_NumberOfPixelsCounted;
}
std::cout<<"TEST2!"<<std::endl;
if( parameters[0] > 0.087)
measure = measure*(parameters[0]*1000);
return measure;
}
/*
* Get both the match Measure and theDerivative Measure
*/
template <class TFixedImage, class TMovingImage>
void
myReg<TFixedImage,TMovingImage>
::GetValueAndDerivative(const TransformParametersType & parameters,
MeasureType & value, DerivativeType & derivative) const
{
itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");
if( !m_GradientImage )
{
itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");
}
FixedImageConstPointer fixedImage = this->GetFixedImage();
if( !fixedImage )
{
itkExceptionMacro( << "Fixed image has not been assigned" );
}
const unsigned int ImageDimension = FixedImageType::ImageDimension;
typedef itk::ImageRegionConstIteratorWithIndex<
FixedImageType> FixedIteratorType;
typedef itk::ImageRegionConstIteratorWithIndex<
ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;
FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );
typename FixedImageType::IndexType index;
MeasureType measure = NumericTraits< MeasureType >::Zero;
m_NumberOfPixelsCounted = 0;
this->SetTransformParameters( parameters );
const unsigned int ParametersDimension = this->GetNumberOfParameters();
derivative = DerivativeType( ParametersDimension );
derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );
ti.GoToBegin();
while(!ti.IsAtEnd())
{
index = ti.GetIndex();
typename Superclass::InputPointType inputPoint;
fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
if( m_FixedImageMask && !m_FixedImageMask->IsInside( inputPoint ) )
{
++ti;
continue;
}
typename Superclass::OutputPointType transformedPoint = m_Transform->TransformPoint( inputPoint );
if( m_MovingImageMask && !m_MovingImageMask->IsInside( transformedPoint ) )
{
++ti;
continue;
}
if( m_Interpolator->IsInsideBuffer( transformedPoint ) )
{
const RealType movingValue = m_Interpolator->Evaluate( transformedPoint );
const TransformJacobianType & jacobian =
m_Transform->GetJacobian( inputPoint );
const RealType fixedValue = ti.Value();
m_NumberOfPixelsCounted++;
const RealType diff = movingValue - fixedValue;
measure += diff * diff;
// Get the gradient by NearestNeighboorInterpolation:
// which is equivalent to round up the point components.
typedef typename Superclass::OutputPointType OutputPointType;
typedef typename OutputPointType::CoordRepType CoordRepType;
typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
MovingImageContinuousIndexType;
MovingImageContinuousIndexType tempIndex;
m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );
typename MovingImageType::IndexType mappedIndex;
for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ )
{
mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
}
const GradientPixelType gradient =
m_GradientImage->GetPixel( mappedIndex );
for(unsigned int par=0; par<ParametersDimension; par++)
{
RealType sum = NumericTraits< RealType >::Zero;
for(unsigned int dim=0; dim<ImageDimension; dim++)
{
sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
}
derivative[par] += sum;
}
}
++ti;
}
if( !m_NumberOfPixelsCounted )
{
itkExceptionMacro(<<"All the points mapped to outside of the moving image");
}
else
{
for(unsigned int i=0; i<ParametersDimension; i++)
{
derivative[i] /= m_NumberOfPixelsCounted;
}
measure /= m_NumberOfPixelsCounted;
}
std::cout<<"TEST 2!"<<std::endl;
if( parameters[0] > 0.087)
measure = measure*(parameters[0]*1000);
value = measure;
}
} // end namespace itk
#endif
"Miller, James V (Research)" <millerjv at crd.ge.com> wrote:
We discussed this a while back. You are essentially asking for a constrained optimization
min F(Image1, Image2, Transform)
s.t. G(Transform) = 0
All the image registration in ITK are unconstrained.
I think there are two ways to do what you want but you will have to write some code.
The first way is to modify or create new optimizers. Most the ITK optimizers are wrappers around vnl optimizers which are wrappers around netlib/lapack type routines. We'd need
to wrap a contrained optimization optimizer and develop an API for specifying constraints. This is probably the best solution architectually.
The second option is to create a new image metric that incorporates your contraint. Here, the metric would incorporate a penalty or additional "cost" for violating the constraints. The metric would evaluate something like:
metric = F(Image1, FImage2, Transfrom) + lambda * G(Transform)
where G(Transform) would be zero for angles below a specified threshold and then increase (linearly?) as the angle exceeds the specified threshold.
The second way requires less code but is architectually fragile.
For example, let's take the MeanSquaresImageToImageMetric. The method GetValue()
takes a set of transformation parameters (vector of parameters) and computes the mean square difference between the images where one of the images is transformed using the specified transformation parameters. Unfortunately, the metric knows nothing about the
type of transformation being used. The metric simply takes the vector of transformation
parameters, and via virtual function call to a transform assigned at runtime, it passes the transformation parameter vector to the transform. (This all works in the registration framework because the registration method will call GetParameters() on the transform and pass them to the optimizer/metric. So the parameter vector is implictly consistent between the transform, metric, and optimizer.)
You could subclass the MeanSquaresImageToImageMetric and add in a penalty term based on the parameter vector passed into routines like GetValue(), GetValueAndDerivative(), etc. But here you would have to peek into specific positions in the parameter vector to check the value of the rotation. This requires knowlege of the type of transform being used. This is why this solution is fragile.
You could make this less fragile by only allowing your new metric to use a certain type of transform. You could override the SetTransform() method in ImageToImageMetric and check that the specified transform is the type of transform needed. Here, you would dynamic cast the pointer to the specified transform to a the specific type of transform your metric uses (for instance Rigid3DTransform, etc.). If the dynamic_cast returns 0, then the wrong type of transform was specified and you would issue an exception.
You could make this less fragile by "deducing" the rotation encapsulated in the transform by transforming the vertices of a simplex (in the GetValue() routines) and groking the rotation applied.
Jim
-----Original Message-----
From: Robert Maroon [mailto:robertmaroon at yahoo.com]
Sent: Friday, October 01, 2004 4:32 PM
To: insight-users at itk.org
Subject: [Insight-users] How to limit rotations when registering in ITK?
Hi all,
I am trying to use the ImageRegistration5 example found in the ITK examples to register a pair of images but I need to restrict the amount of rotation allowed. I have looked though the examples and documentation and I can't seem to find a way to constrain the angles (say less than x degrees) in which the registration can look. I has looking particular to see if there is anything besides the scaling and MaximumStepLengths that can be set for the optimizer that might let me acheive this.
Thanks!
Robert
---------------------------------
Do you Yahoo!?
Yahoo! Mail - You care about security. So do we.
---------------------------------
Do you Yahoo!?
vote.yahoo.com - Register online to vote today!_______________________________________________
Insight-users mailing list
Insight-users at itk.org
http://www.itk.org/mailman/listinfo/insight-users
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20041008/6deafa55/attachment-0001.htm
More information about the Insight-users
mailing list