[Insight-users] application of ModelToImgBasedReg. : PCA to Image Registration

Luis Ibanez luis . ibanez at kitware . com
Sun, 05 Oct 2003 20:08:16 -0400


This is a multi-part message in MIME format.
--------------020805060106080509080009
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit


Hi Steven,

If you have a PCA model of the cartilage, the set of paramaters for the
registartion will be a set of K values:

                                       K = M  + N

where:
 
N is the number of principal components you
    decide to retain from the PCA analysis
     (I assume that's what you meant when you say that
      you need some "extra variance parameter(s) to be
     perturbed")

and

M  is the number of parameters of the geometrical
      transformation that places your PCA model in space.

Let's say for example that you decided to retain the  five principal 
components
of the PCA analysis, and that you use a rigid registration in 3D.

A rigid transform will require 6 parametes ( 3 for translation and 3 for 
rotation).

Your total set of parameters will then be:  11 = 5(from PCA) + 6(from 
transform).

Then, you have to be very careful with the definition of the Metric that 
will
evaluate the matching of your PCA model to the Image.

By default, current ITK registration metrics assume that the parameters are
those of the transform.  In your case however, you have to take the 
vector of
parameters and separate some that will go to the PCA model  (5 of the 
parameters)
and some that will go to the rigid transform (6 -> for the 6 degrees of 
freedom).

Is your PCA model formed with points in the surface of the cartilage ?
If so, you should register your model agains the GradientMagnitude of the
input image, instead of the original image.  

You could use as a metric the simple sum of Image values (in this case
Gradient Magnitude values) on the positions indicated by the points of the
PCA model. In this way, when the points of the PCA model (for a particular
set of modes values) is perfectly superimposed to the edge of the cartilage,
the sum of values from the Gradient Magnitude will be maximized.

Please look at the attached file that implements (roughtly) the GetValue()
method of a PCA to Image metric.

It assumes that your have a C++ class that represents the PCA Model, and
this class responds to the interface:

  Begin(), End() : returning iterators to the list of points in the PCA 
model
  (you could use an itk::PointSet for the internal representation,... or an
   itk::SpatialObject).

  SetPCAModes( std::vector<double> )  where you can set the current values
    of the modes associated with the principal components.

This code will give you an idea of how to implement your PCA to Image
metric.

Once you are done with the metric you can use all the other components
of the registration framework and configure a full registration method.

It souldn't be too hard to do multiresolution also, by modifying the
MultiRes Image Registration method for accepting one PCAModel
class and an Image Pyramid. (there may be some Coffee involved though...)



Please let us know if you have further questions or experience
difficulties implementing the PCA to image Metric.


Regards,


   Luis


--------------------------------------
sdor5151 at mail . usyd . edu . au wrote:

>Dear ITK users,
>
>  I've been looking at methods of segmenting the knee cartilage with Active
>Shape Models.
>
>  I have built a statistical shape model of the cartilage from a training set.
>
>  Now I was wondering whether ModelToImageBasedRegistration.cxx can be used to
>actually match the model to a target MRI scan. (since there isn't such thing as
>a Multiresolution search for a model in an image). 
>
>  If so, how can I modify the ParametersType in GetValue() so that there is an
>extra variance parameter to be perturbed when doing the search?
>
>  Thankyou,
>
>Steven
>
>
>
>-------------------------------------------------
>This mail sent through IMP: www-mail.usyd.edu.au
>_______________________________________________
>Insight-users mailing list
>Insight-users at itk . org
>http://www . itk . org/mailman/listinfo/insight-users
>
>  
>
----------------------------------------------------------------------




--------------020805060106080509080009
Content-Type: text/plain;
 name="itkPCAToImageMetric.txx"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="itkPCAToImageMetric.txx"


template <class TFixedPCAModel, class TMovingImage> 
typename PCAModelToImageMetric<TFixedPCAModel,TMovingImage>::MeasureType
PCAModelToImageMetric<TFixedPCAModel,TMovingImage>
::GetValue( const TransformParametersType & parameters ) const
{

  itkDebugMacro("GetValue( " << parameters << " ) ");

  FixedPCAModelConstPointer fixedPCAModel = this->GetFixedPCAModel();

  if( !fixedPCAModel ) 
    {
    itkExceptionMacro( << "Fixed PCAModel has not been assigned" );
    }

  // HERE Take the parameters that control the modes of the PCA model
  std::vector<double> PCAmodes;
  for(unsigned int i=0; i<NumberOfPCAModes; i++)
    {
    PCAModes.push_back( parameters[i] );
    }

  fixedPCAModel->SetModes( PCAmodes );


  FixedPCAModelIteratorType PCAPoint = fixedPCAModel->Begin();
  FixedPCAModelIteratorType PCALast  = fixedPCAModel->End();


  MeasureType measure = NumericTraits< MeasureType >::Zero;

  m_NumberOfPixelsCounted = 0;

  // HERE copy the remaining parameters for the transform
  for(unsigned int i=0; i<NumberOFTransformParamters; i++)
    {
    transformParameters[i] =  parameters[i+NumberOfPCAModes] );
    }


  this->SetTransformParameters( transformParameters );

  while(!ti.IsAtEnd())
    {

    typename Superclass::InputPointType inputPoint = *PCAPoint;

    typename Superclass::OutputPointType transformedPoint = m_Transform->TransformPoint( inputPoint );

    if( m_Interpolator->IsInsideBuffer( transformedPoint ) )
      {
      const RealType movingValue  = m_Interpolator->Evaluate( transformedPoint );
      m_NumberOfPixelsCounted++;
      measure += movingValue;
      }
    ++PCAPoint;
    }

  if( !m_NumberOfPixelsCounted )
    {
    itkExceptionMacro(<<"All the points mapped to outside of the moving image");
    }

  return measure;

}





--------------020805060106080509080009--