<DIV>Hello Robert,</DIV>
<DIV> </DIV>
<DIV>The LBFGS<STRONG>B</STRONG>Optimizer 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.</DIV>
<DIV> </DIV>
<DIV>Edi</DIV>
<DIV> </DIV>
<DIV> </DIV>
<DIV><BR><BR><B><I>Robert Maroon <robertmaroon@yahoo.com></I></B> wrote:</DIV>
<BLOCKQUOTE class=replbq style="PADDING-LEFT: 5px; MARGIN-LEFT: 5px; BORDER-LEFT: #1010ff 2px solid">
<DIV><FONT face="times new roman" size=3>Thanks for the advice James,</FONT></DIV>
<DIV><FONT face="times new roman" size=3></FONT> </DIV>
<DIV><FONT face="times new roman" size=3>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</FONT><FONT face=Verdana> 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. </FONT></DIV>
<DIV><FONT face=Verdana></FONT> </DIV>
<DIV><FONT face=Verdana>Thanks in advance! I would really appreciate some help!</FONT></DIV>
<DIV><FONT face=Verdana></FONT> </DIV>
<DIV><FONT face=Verdana>In my main project code I call on my class as follows:</FONT></DIV>
<DIV><FONT face=Verdana></FONT> </DIV>
<DIV><FONT face=Verdana>I include:</FONT></DIV>
<DIV>#include "myReg.h"</DIV>
<DIV><FONT face=Verdana></FONT> </DIV>
<DIV><FONT face=Verdana>And then to declare the metric I do:</FONT></DIV>
<DIV>typedef itk::myReg<FixedImageType,MovingImageType> MetricType;</DIV>
<DIV> </DIV>
<DIV>The rest of the file is basically ImageRegistration5.cxx</DIV>
<DIV><FONT face=Verdana></FONT> </DIV>
<DIV><FONT face=Verdana></FONT> </DIV>
<DIV>#ifndef __myReg_h<BR>#define __myReg_h</DIV>
<DIV>#include "itkImageToImageMetric.h"<BR>#include "itkCovariantVector.h"<BR>#include "itkMeanSquaresImageToImageMetric.h"<BR>#include "itkPoint.h"</DIV>
<DIV><BR>namespace itk<BR>{<BR><BR>template < class TFixedImage, class TMovingImage > <BR>class ITK_EXPORT myReg : <BR> public MeanSquaresImageToImageMetric< TFixedImage, TMovingImage><BR>{</DIV>
<DIV> /** Get the value for single valued optimizers. */<BR> MeasureType GetValue( const TransformParametersType & parameters ) const;</DIV>
<DIV> /** Get value and derivatives for multiple valued optimizers. */<BR> void GetValueAndDerivative( const TransformParametersType & parameters,<BR> MeasureType& Value, DerivativeType& Derivative ) const;</DIV>
<DIV>protected:<BR> myReg();<BR> virtual ~myReg() {};</DIV>
<DIV>};</DIV>
<DIV>} // end namespace itk</DIV>
<DIV><BR>#endif</DIV>
<DIV><FONT face=Verdana></FONT> </DIV>
<DIV><FONT face=Verdana></FONT> </DIV>
<DIV><FONT face=Verdana>Here's the TXX file:<BR><BR>#ifndef _myReg_txx<BR>#define _myReg_txx</FONT></DIV>
<DIV><FONT face=Verdana>#include "myReg.h"<BR>#include "itkImageRegionConstIteratorWithIndex.h"</FONT></DIV>
<DIV><FONT face=Verdana>namespace itk<BR>{</FONT></DIV>
<DIV><FONT face=Verdana>/*<BR> * Constructor<BR> */<BR>template <class TFixedImage, class TMovingImage> <BR>myReg<TFixedImage,TMovingImage><BR>::myReg()<BR>{<BR> std::cout<<"TEST 1!"<<std::endl;<BR> itkDebugMacro("Constructor");<BR>}</FONT></DIV>
<DIV><FONT face=Verdana>/*<BR> * Get the match Measure<BR> */<BR>template <class TFixedImage, class TMovingImage> <BR>typename myReg<TFixedImage,TMovingImage>::MeasureType<BR>myReg<TFixedImage,TMovingImage><BR>::GetValue( const TransformParametersType & parameters ) const<BR>{</FONT></DIV>
<DIV><FONT face=Verdana> itkDebugMacro("GetValue( " << parameters << " ) ");</FONT></DIV>
<DIV><FONT face=Verdana> FixedImageConstPointer fixedImage = this->GetFixedImage();</FONT></DIV>
<DIV><FONT face=Verdana> if( !fixedImage ) <BR> {<BR> itkExceptionMacro( << "Fixed image has not been assigned" );<BR> }</FONT></DIV>
<DIV><FONT face=Verdana> typedef itk::ImageRegionConstIteratorWithIndex<FixedImageType> FixedIteratorType;</FONT></DIV><FONT face=Verdana>
<DIV><BR> FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );</DIV>
<DIV> typename FixedImageType::IndexType index;</DIV>
<DIV> MeasureType measure = NumericTraits< MeasureType >::Zero;</DIV>
<DIV> m_NumberOfPixelsCounted = 0;</DIV>
<DIV> this->SetTransformParameters( parameters );</DIV>
<DIV> while(!ti.IsAtEnd())<BR> {</DIV>
<DIV> index = ti.GetIndex();<BR> <BR> typename Superclass::InputPointType inputPoint;<BR> fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );</DIV>
<DIV> if( m_FixedImageMask && !m_FixedImageMask->IsInside( inputPoint ) )<BR> {<BR> ++ti;<BR> continue;<BR> }</DIV>
<DIV> typename Superclass::OutputPointType transformedPoint = m_Transform->TransformPoint( inputPoint );</DIV>
<DIV> if( m_MovingImageMask && !m_MovingImageMask->IsInside( transformedPoint ) )<BR> {<BR> ++ti;<BR> continue;<BR> }</DIV>
<DIV> if( m_Interpolator->IsInsideBuffer( transformedPoint ) )<BR> {<BR> const RealType movingValue = m_Interpolator->Evaluate( transformedPoint );<BR> const RealType fixedValue = ti.Get();<BR> m_NumberOfPixelsCounted++;<BR> const RealType diff = movingValue - fixedValue; <BR> measure += diff * diff; <BR> }</DIV>
<DIV> ++ti;<BR> }</DIV>
<DIV> if( !m_NumberOfPixelsCounted )<BR> {<BR> itkExceptionMacro(<<"All the points mapped to outside of the moving image");<BR> }<BR> else<BR> {<BR> measure /= m_NumberOfPixelsCounted;<BR> }</DIV>
<DIV> std::cout<<"TEST2!"<<std::endl;</DIV>
<DIV> if( parameters[0] > 0.087)<BR> measure = measure*(parameters[0]*1000);</DIV>
<DIV> return measure;</DIV>
<DIV>}</DIV>
<DIV> </DIV>
<DIV>/*<BR> * Get both the match Measure and theDerivative Measure <BR> */<BR>template <class TFixedImage, class TMovingImage> <BR>void<BR>myReg<TFixedImage,TMovingImage><BR>::GetValueAndDerivative(const TransformParametersType & parameters, <BR> MeasureType & value, DerivativeType & derivative) const<BR>{</DIV>
<DIV> itkDebugMacro("GetValueAndDerivative( " << parameters << " ) ");</DIV>
<DIV> if( !m_GradientImage )<BR> {<BR> itkExceptionMacro(<<"The gradient image is null, maybe you forgot to call Initialize()");<BR> }</DIV>
<DIV> FixedImageConstPointer fixedImage = this->GetFixedImage();</DIV>
<DIV> if( !fixedImage ) <BR> {<BR> itkExceptionMacro( << "Fixed image has not been assigned" );<BR> }</DIV>
<DIV> const unsigned int ImageDimension = FixedImageType::ImageDimension;</DIV>
<DIV> typedef itk::ImageRegionConstIteratorWithIndex<<BR> FixedImageType> FixedIteratorType;</DIV>
<DIV> typedef itk::ImageRegionConstIteratorWithIndex<<BR> ITK_TYPENAME Superclass::GradientImageType> GradientIteratorType;</DIV>
<DIV><BR> FixedIteratorType ti( fixedImage, this->GetFixedImageRegion() );</DIV>
<DIV> typename FixedImageType::IndexType index;</DIV>
<DIV> MeasureType measure = NumericTraits< MeasureType >::Zero;</DIV>
<DIV> m_NumberOfPixelsCounted = 0;</DIV>
<DIV> this->SetTransformParameters( parameters );</DIV>
<DIV> const unsigned int ParametersDimension = this->GetNumberOfParameters();<BR> derivative = DerivativeType( ParametersDimension );<BR> derivative.Fill( NumericTraits<ITK_TYPENAME DerivativeType::ValueType>::Zero );</DIV>
<DIV> ti.GoToBegin();</DIV>
<DIV> while(!ti.IsAtEnd())<BR> {</DIV>
<DIV> index = ti.GetIndex();<BR> <BR> typename Superclass::InputPointType inputPoint;<BR> fixedImage->TransformIndexToPhysicalPoint( index, inputPoint );</DIV>
<DIV> if( m_FixedImageMask && !m_FixedImageMask->IsInside( inputPoint ) )<BR> {<BR> ++ti;<BR> continue;<BR> }</DIV>
<DIV> typename Superclass::OutputPointType transformedPoint = m_Transform->TransformPoint( inputPoint );</DIV>
<DIV> if( m_MovingImageMask && !m_MovingImageMask->IsInside( transformedPoint ) )<BR> {<BR> ++ti;<BR> continue;<BR> }</DIV>
<DIV> if( m_Interpolator->IsInsideBuffer( transformedPoint ) )<BR> {<BR> const RealType movingValue = m_Interpolator->Evaluate( transformedPoint );</DIV>
<DIV> const TransformJacobianType & jacobian =<BR> m_Transform->GetJacobian( inputPoint ); </DIV>
<DIV> <BR> const RealType fixedValue = ti.Value();<BR> m_NumberOfPixelsCounted++;</DIV>
<DIV> const RealType diff = movingValue - fixedValue; <BR> <BR> measure += diff * diff;</DIV>
<DIV> // Get the gradient by NearestNeighboorInterpolation: <BR> // which is equivalent to round up the point components.<BR> typedef typename Superclass::OutputPointType OutputPointType;<BR> typedef typename OutputPointType::CoordRepType CoordRepType;<BR> typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension><BR> MovingImageContinuousIndexType;</DIV>
<DIV> MovingImageContinuousIndexType tempIndex;<BR> m_MovingImage->TransformPhysicalPointToContinuousIndex( transformedPoint, tempIndex );</DIV>
<DIV> typename MovingImageType::IndexType mappedIndex; <BR> for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ )<BR> {<BR> mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );<BR> }</DIV>
<DIV> const GradientPixelType gradient = <BR> m_GradientImage->GetPixel( mappedIndex );</DIV>
<DIV> for(unsigned int par=0; par<ParametersDimension; par++)<BR> {<BR> RealType sum = NumericTraits< RealType >::Zero;<BR> for(unsigned int dim=0; dim<ImageDimension; dim++)<BR> {<BR> sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];<BR> }<BR> derivative[par] += sum;<BR> }<BR> }</DIV>
<DIV> ++ti;<BR> }</DIV>
<DIV> if( !m_NumberOfPixelsCounted )<BR> {<BR> itkExceptionMacro(<<"All the points mapped to outside of the moving image");<BR> }<BR> else<BR> {<BR> for(unsigned int i=0; i<ParametersDimension; i++)<BR> {<BR> derivative[i] /= m_NumberOfPixelsCounted;<BR> }<BR> measure /= m_NumberOfPixelsCounted;<BR> }</DIV>
<DIV> std::cout<<"TEST 2!"<<std::endl;</DIV>
<DIV> if( parameters[0] > 0.087)<BR> measure = measure*(parameters[0]*1000);</DIV>
<DIV> value = measure;</DIV>
<DIV>}</DIV>
<DIV>} // end namespace itk</DIV>
<DIV><BR>#endif<BR></DIV>
<DIV> </DIV>
<DIV> </DIV></FONT>
<DIV><FONT face=Verdana></FONT> </DIV>
<DIV><FONT face=Verdana></FONT> </DIV>
<P><FONT face=Verdana></FONT> </P>
<P><FONT face=Verdana></FONT> </P>
<P><FONT face=Verdana></FONT> </P>
<DIV><BR><BR><B><I>"Miller, James V (Research)" <millerjv@crd.ge.com></I></B> wrote:</DIV>
<BLOCKQUOTE class=replbq style="PADDING-LEFT: 5px; MARGIN-LEFT: 5px; BORDER-LEFT: #1010ff 2px solid">
<META content="MSHTML 6.00.2800.1276" name=GENERATOR>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>We discussed this a while back. You are essentially asking for a constrained optimization</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=677103113-04102004> <FONT face=Verdana color=#0000ff size=2>min F(Image1, Image2, Transform)</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004></SPAN><SPAN class=677103113-04102004> <FONT face=Verdana color=#0000ff size=2>s.t. G(Transform) = 0</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>All the image registration in ITK are unconstrained.</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>I think there are two ways to do what you want but you will have to write some code.</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>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 </FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>to wrap a contrained optimization optimizer and develop an API for specifying constraints. This is probably the best solution architectually.</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>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:</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=677103113-04102004> <FONT face=Verdana color=#0000ff size=2>metric = F(Image1, FImage2, Transfrom) + lambda * G(Transform)</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>where G(Transform) would be zero for angles below a specified threshold and then increase (linearly?) as the angle exceeds the specified threshold.</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>The second way requires less code but is architectually fragile. </FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>For example, let's take the MeanSquaresImageToImageMetric. The method GetValue()</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>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 </FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>type of transformation being used. The metric simply takes the vector of transformation </FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>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.)</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>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.</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>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.</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>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.</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN> </DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2>Jim</FONT></SPAN></DIV>
<DIV><SPAN class=677103113-04102004><FONT face=Verdana color=#0000ff size=2></FONT></SPAN> </DIV>
<BLOCKQUOTE>
<DIV class=OutlookMessageHeader dir=ltr align=left><FONT face=Tahoma size=2>-----Original Message-----<BR><B>From:</B> Robert Maroon [mailto:robertmaroon@yahoo.com]<BR><B>Sent:</B> Friday, October 01, 2004 4:32 PM<BR><B>To:</B> insight-users@itk.org<BR><B>Subject:</B> [Insight-users] How to limit rotations when registering in ITK?<BR><BR></FONT></DIV>
<DIV>Hi all,</DIV>
<DIV> </DIV>
<DIV>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.</DIV>
<DIV> </DIV>
<DIV>Thanks!<BR></DIV>
<DIV>Robert</DIV>
<P>
<HR SIZE=1>
Do you Yahoo!?<BR><A href="http://us.rd.yahoo.com/mail_us/taglines/security/*http://promotions.yahoo.com/new_mail/static/protection.html">Yahoo! Mail</A> - You care about security. So do we.</BLOCKQUOTE></BLOCKQUOTE>
<P>
<HR SIZE=1>
Do you Yahoo!?<BR><A href="http://vote.yahoo.com/">vote.yahoo.com</A> - Register online to vote today!_______________________________________________<BR>Insight-users mailing list<BR>Insight-users@itk.org<BR>http://www.itk.org/mailman/listinfo/insight-users<BR></BLOCKQUOTE>