itkPathFunctions.h

Go to the documentation of this file.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkPathFunctions.h,v $
  Language:  C++
  Date:      $Date: 2009-04-06 00:08:10 $
  Version:   $Revision: 1.9 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/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 __itkPathFunctions_h
#define __itkPathFunctions_h


#include "itkPath.h"
#include "itkChainCodePath.h"
#include "itkFourierSeriesPath.h"
#include "itkOffset.h"
#include <math.h>

namespace itk
{
template <class TChainCodePath, class TPathInput>
void MakeChainCodeTracePath( TChainCodePath & chainPath,
                             const TPathInput & inPath,
                             bool restrictMovement = false )
{
  typedef typename TChainCodePath::OffsetType OffsetType;
  typedef typename TChainCodePath::InputType  ChainInputType;
  typedef typename TChainCodePath::OutputType ChainOutputType;
  typedef typename TPathInput::InputType      InPathInputType;
  typedef typename TPathInput::OutputType     InPathOutputType;

  OffsetType offset, tempOffset, zeroOffset;
  InPathInputType inPathInput;
  int dimension = OffsetType::GetOffsetDimension();
  
  zeroOffset.Fill(0);
  
  chainPath.Clear();
  inPathInput = inPath.StartOfInput();
  chainPath.SetStart(  inPath.EvaluateToIndex( inPathInput )  );
  
  for(ChainInputType chainInput=0;;)
    {
    offset  = inPath.IncrementInput(inPathInput);
    if( zeroOffset == offset ) { break; }
    
    if( ! restrictMovement )
      {
      chainPath.InsertStep( chainInput++, offset );
      }
    else
      {
      for( int d=0; d<dimension; d++ )
        {
        tempOffset.Fill(0);
        tempOffset[d] = offset[d];
        chainPath.InsertStep( chainInput++, tempOffset );
        }
      }
    }
}

template <class TFourierSeriesPath, class TChainCodePath>
void MakeFourierSeriesPathTraceChainCode( TFourierSeriesPath & FSPath,
                             const TChainCodePath & chainPath,
                             unsigned int numHarmonics = 8 )
{
  typedef typename TFourierSeriesPath::IndexType  IndexType;
  typedef typename TFourierSeriesPath::OffsetType OffsetType;
  typedef typename TFourierSeriesPath::VectorType VectorType;

  typedef typename TFourierSeriesPath::InputType  FSInputType;
  typedef typename TFourierSeriesPath::OutputType FSOutputType;
  typedef typename TChainCodePath::InputType      ChainInputType;
  typedef typename TChainCodePath::OutputType     ChainOutputType;
  
  IndexType   index;
  VectorType  indexVector;
  VectorType  cosCoefficient;
  VectorType  sinCoefficient;
  FSInputType theta;
  int         dimension =     OffsetType::GetOffsetDimension();
  unsigned    numSteps  =     chainPath.NumberOfSteps();
  
  const double PI = 4.0 * vcl_atan( 1.0 );

  FSPath.Clear();
  
  // Adjust our private copy of numHarmonics if necessary
  if( numHarmonics <= 1 )
    numHarmonics = 2;
  else if( numHarmonics*2 > numSteps )
    numHarmonics = numSteps / 2;
  
  for( unsigned n=0; n<numHarmonics; n++ )
    {
    index = chainPath.GetStart();
    cosCoefficient.Fill(0.0);
    sinCoefficient.Fill(0.0);
    
    for( ChainInputType step=0; step<numSteps; step++ )
      {
      index += chainPath.Evaluate( step );
      theta = 2 * n * PI * ( double(step+1)) / numSteps;
      
      // turn the current index into a vector
      for( int d=0; d<dimension; d++ )
        {
        indexVector[d] = index[d];
        }
      cosCoefficient += indexVector * (vcl_cos(theta)/numSteps);
      sinCoefficient += indexVector * (vcl_sin(theta)/numSteps);
      }
    
    FSPath.AddHarmonic( cosCoefficient, sinCoefficient );
    }
}

} // end namespace itk

#endif