freDemonRegistrationProcessor.txx

Go to the documentation of this file.

/*=========================================================================
 
  Program:   F.R.E.E. - flexible registration evaluation engine
  Version:   v.1.0.0
  Date:      $Date: 2006/09/01 12:00:00 $
  Module:    $RCSfile: freDemonRegistrationProcessor.h,v $
  Language:  C++
 
 
 
  Copyright (c) 2007 Ralf o Floca (Department of Medical Informatics,
  Institute for Medical Biometry and Informatics, University of Heidelberg,
  Germany). All rights reserved.
  See FREECopyright.txt or http://www.mi.med.uni-hd.de/free/copyright.htm
  for details.
 
     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.
 
=========================================================================*/
#ifndef __freDemonRegistrationProcessor_txx
#define __freDemonRegistrationProcessor_txx

#include "freDemonRegistrationProcessor.h"
#include "freOptimizerControllerBase.h"
#include "freTransformSetupAdaptor.h"

#include "itkPasteImageFilter.h"

namespace FREE
{



template <unsigned int VImageDimension>
void
DemonRegistrationProcessor<VImageDimension>::
ComputeRegistration()
{
    long optimizerObserverID = 0;
    try
    {
        Superclass::ComputeRegistration();

        typename HistoMatchFilterType::Pointer histoMatch = HistoMatchFilterType::New();

        ComponentSetup* pSetup = this->GetComponentSetup();
        int iNrOfIterations;
        if(!pSetup->Parameters().GetParameterValue(cParam_Iterations,iNrOfIterations))
            throwExceptionMacro("Missing parameter: " << cParam_Iterations);

        double dStdDev;
        if(!pSetup->Parameters().GetParameterValue(cParam_StdDev,dStdDev))
            throwExceptionMacro("Missing parameter: " << cParam_StdDev);

        bool bHistoMatch;
        if(!pSetup->Parameters().GetParameterValue(cParam_HistoMatch,bHistoMatch))
            throwExceptionMacro("Missing parameter: " << cParam_HistoMatch);

        long lHistoLevels;
        if(!pSetup->Parameters().GetParameterValue(cParam_HistoLevels,lHistoLevels))
            throwExceptionMacro("Missing parameter: " << cParam_HistoLevels);

        long lHistoPoints;
        if(!pSetup->Parameters().GetParameterValue(cParam_HistoMatchPoints,lHistoPoints))
            throwExceptionMacro("Missing parameter: " << cParam_HistoMatchPoints);

        bool bHistoThreshold;
        if(!pSetup->Parameters().GetParameterValue(cParam_HistoThreshold,bHistoThreshold))
            throwExceptionMacro("Missing parameter: " << cParam_HistoThreshold);

        //Extract needed ImageParts;
        typename ImageType::Pointer smpMovingPart = this->m_smpMovingImage;
        typename ImageType::Pointer smpFixedPart = this->m_smpFixedImage;

        //Connect images
        if(bHistoMatch)
        {
            if(this->fnOnProgress.IsNotNull())
                this->fnOnProgress->Execute(RSTInitProcessor, "Histogram match enabled",this);

            histoMatch->SetInput( smpMovingPart );
            histoMatch->SetReferenceImage( smpFixedPart );

            histoMatch->SetNumberOfHistogramLevels( lHistoLevels );
            histoMatch->SetNumberOfMatchPoints( lHistoPoints );
            if ( bHistoThreshold )
                histoMatch->ThresholdAtMeanIntensityOn();

            //Calculate histogram matching
            histoMatch->Update();
            //Connect it to registration
            smpRegistration->SetFixedImage(smpFixedPart);
            smpRegistration->SetMovingImage(histoMatch->GetOutput());
        }
        else
        {
            smpRegistration->SetFixedImage(smpFixedPart);
            smpRegistration->SetMovingImage(smpMovingPart);
        };

        smpRegistration->SetNumberOfIterations( iNrOfIterations );
        smpRegistration->SetStandardDeviations( dStdDev );

        //Register observer
        m_IterationObserver->fnOnNotify = m_IterationNotificationEvent;
        optimizerObserverID = this->smpRegistration->AddObserver(itk::IterationEvent(), this->m_IterationObserver);
    }
    catchAllNPassMacro("Unknown Error while preparing the registration.");

    try
    {
        //Call the OnNextLevel one time before starting the registration.
        //even a singel resolution Registration has at least one level to notify
        if(this->fnOnNextLevel.IsNotNull())
            this->fnOnNextLevel->Execute(0,this);

        this->smpRegistration->Update();

        this->smpRegistration->RemoveObserver(optimizerObserverID);
    }
    catchAllNPassMacro("Unknown Error while itk processing the registration.");
};

template <unsigned int VImageDimension>
DemonRegistrationProcessor<VImageDimension>::
DemonRegistrationProcessor()
{
    smpRegistration = itk::DemonsRegistrationFilter< ImageType, ImageType, TransformFieldType >::New();
    smpInterpolate = NULL;
    m_IterationObserver = IterationObserver::New();
    m_IterationNotificationEvent = NotificationEvent<DemonRegistrationProcessor>::New(this, &Self::OnNewIteration);
};

template <unsigned int VImageDimension>
DemonRegistrationProcessor<VImageDimension>::
~DemonRegistrationProcessor()
{}
;

template <unsigned int VImageDimension>
typename DemonRegistrationProcessor<VImageDimension>::TransformFieldPointer
DemonRegistrationProcessor<VImageDimension>::
ComputeTransformationField(const PointType& origin, const RegionType& region, const SpacingType& spacing)
{
    if(this->fnOnProgress.IsNotNull())
        this->fnOnProgress->Execute(RSTFinalizing, "Computing transformation field",this);

    TransformFieldPointer smpDemField = smpRegistration->GetDeformationField();

    if (!( (region==smpDemField->GetLargestPossibleRegion())&
           (spacing==smpDemField->GetSpacing())&
           (origin==smpDemField->GetOrigin()) ) )
    { // have to resize the field.

        //Create the field in correct size
        typename TransformFieldType::Pointer smpField = TransformFieldType::New();
        smpField->SetOrigin(origin);
        smpField->SetSpacing(spacing);
        smpField->SetRegions(region);
        smpField->Allocate();
        //Reset the transformation field to no transformation (all vectors components are 0)
        typename TransformFieldType::PixelType nullVector;
        for (int iIndex=0; iIndex < VImageDimension; iIndex++)
            nullVector[iIndex]=0;

        // iterator for the output image
        itk::ImageRegionIteratorWithIndex<TransformFieldType> outputIt(
          smpField, smpField->GetLargestPossibleRegion() );

        // interpolator
        typedef itk::LinearInterpolateImageFunction<TransformFieldType,ScalarType> InterpolatorType;
        typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
        interpolator->SetInputImage(smpDemField);

        typename ImageType::IndexType index;
        typename ImageType::PointType point;
        typename TransformFieldType::PixelType value;

        while( !outputIt.IsAtEnd() )
        {
          // get the output image index
          index = outputIt.GetIndex();
          smpField->TransformIndexToPhysicalPoint( index, point );

          if( interpolator->IsInsideBuffer( point ) )
          {
            value = interpolator->Evaluate( point );
            outputIt.Set( value );
          }
          else
          {
            outputIt.Set( nullVector );
          }  

          ++outputIt;
        }

        return smpField;
    }
    else
        return smpDemField;
};

template <unsigned int VImageDimension>
void
DemonRegistrationProcessor<VImageDimension>::
OnNewIteration(void* pSender, long threadID)
{
    if (this->fnOnNextIteration)
    {
        //(*fnOnNextIteration)(lIteration, value, parameters,this, 0);
    };
};

}//End of Namespace free

#endif

---