cmtkVoxelMatchingElasticFunctional.h

Go to the documentation of this file.

/*
//
//  Copyright 1997-2010 Torsten Rohlfing
//
//  Copyright 2004-2011 SRI International
//
//  This file is part of the Computational Morphometry Toolkit.
//
//  http://www.nitrc.org/projects/cmtk/
//
//  The Computational Morphometry Toolkit is free software: you can
//  redistribute it and/or modify it under the terms of the GNU General Public
//  License as published by the Free Software Foundation, either version 3 of
//  the License, or (at your option) any later version.
//
//  The Computational Morphometry Toolkit is distributed in the hope that it
//  will be useful, but WITHOUT ANY WARRANTY; without even the implied
//  warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License along
//  with the Computational Morphometry Toolkit.  If not, see
//  <http://www.gnu.org/licenses/>.
//
//  $Revision: 2752 $
//
//  $LastChangedDate: 2011-01-17 11:33:31 -0800 (Mon, 17 Jan 2011) $
//
//  $LastChangedBy: torstenrohlfing $
//
*/

#ifndef __cmtkVoxelMatchingElasticFunctional_h_included_
#define __cmtkVoxelMatchingElasticFunctional_h_included_

#include <cmtkconfig.h>

#include <Registration/cmtkVoxelMatchingFunctional.h>

#include <Base/cmtkVector.h>
#include <Base/cmtkSplineWarpXform.h>
#include <Base/cmtkJointHistogram.h>
#include <Base/cmtkUniformVolume.h>
#include <Base/cmtkMacros.h>
#include <Base/cmtkMathUtil.h>

#include <System/cmtkException.h>

#include <cassert>
#include <vector>

#ifdef HAVE_IEEEFP_H
#  include <ieeefp.h>
#endif

#ifdef HAVE_VALUES_H
#  include <values.h>
#endif

namespace
cmtk
{


class VoxelMatchingElasticFunctional : 
  public VoxelMatchingFunctional 
{
public:
  typedef VoxelMatchingElasticFunctional Self;

  typedef SmartPointer<Self> SmartPtr;

  typedef VoxelMatchingFunctional Superclass;

  virtual void SetWarpXform( SplineWarpXform::SmartPtr& warp ) = 0;

  virtual void SetForceOutside( const bool flag = true, const Types::DataItem value = 0 ) = 0;

  virtual ~VoxelMatchingElasticFunctional ();

protected:
  VoxelMatchingElasticFunctional( UniformVolume::SmartPtr& reference, UniformVolume::SmartPtr& floating );

  cmtkGetSetMacroDefault(bool,AdaptiveFixParameters,true);

  cmtkGetSetMacro(double,AdaptiveFixThreshFactor);

  cmtkGetSetMacroString(ActiveCoordinates);

  cmtkGetSetMacroDefault(double,JacobianConstraintWeight,0);

  cmtkGetSetMacroDefault(double,RigidityConstraintWeight,0);

  cmtkGetSetMacro(DataGrid::SmartPtr,RigidityConstraintMap);

  cmtkGetSetMacroDefault(double,GridEnergyWeight,0);

  cmtkGetSetMacroDefault(bool,Regularize,false);

  bool WarpNeedsFixUpdate;

  JointHistogram<unsigned int>::SmartPtr ConsistencyHistogram;

  size_t Dim;

  std::vector<Types::Coordinate> StepScaleVector;

  DataGrid::RegionType *VolumeOfInfluence;

  Vector3D ReferenceFrom;

  Vector3D ReferenceTo;

  Vector3D *VectorCache;
};

template<class W> 
class VoxelMatchingElasticFunctional_WarpTemplate :
  public VoxelMatchingElasticFunctional 
{
public:
  typedef VoxelMatchingElasticFunctional_WarpTemplate<W> Self;

  typedef SmartPointer<Self> SmartPtr;

  typedef VoxelMatchingElasticFunctional Superclass;

  typename W::SmartPtr Warp;

protected:
  typename W::SmartPtr InverseTransformation;

  double InverseConsistencyWeight;

public:
  void SetInverseTransformation( typename W::SmartPtr& inverseTransformation ) 
  {
    this->InverseTransformation = W::SmartPtr::DynamicCastFrom( inverseTransformation );
  }

  void SetInverseConsistencyWeight( const double inverseConsistencyWeight ) 
  {
    this->InverseConsistencyWeight = inverseConsistencyWeight;
  }

protected:

        typename Self::ReturnType WeightedTotal( const typename Self::ReturnType metric, const W* warp ) const 
  {
    double result = metric;
    if ( this->m_JacobianConstraintWeight > 0 ) 
      {
      result -= this->m_JacobianConstraintWeight * warp->GetJacobianConstraint();
      } 
    
    if ( this->m_RigidityConstraintWeight > 0 ) 
      {
      if ( this->m_RigidityConstraintMap )
        {
        result -= this->m_RigidityConstraintWeight * warp->GetRigidityConstraint( this->m_RigidityConstraintMap );
        }
      else
        {
        result -= this->m_RigidityConstraintWeight * warp->GetRigidityConstraint();
        }
      } 
    
    if ( this->m_GridEnergyWeight > 0 ) 
      {
      result -= this->m_GridEnergyWeight * warp->GetGridEnergy();
      }
    
    if ( !finite( result ) ) 
      return -FLT_MAX;
    
    if ( this->m_MatchedLandmarkList ) 
      {
      result -= this->m_LandmarkErrorWeight * warp->GetLandmarksMSD( this->m_MatchedLandmarkList );
      }

    if ( InverseTransformation ) 
      {
      result -= this->InverseConsistencyWeight * warp->GetInverseConsistencyError( this->InverseTransformation, this->ReferenceGrid );
      }
    
    return static_cast<typename Self::ReturnType>( result );
  }
  
  void WeightedDerivative( double& lower, double& upper, W& warp, const int param, const Types::Coordinate step ) const;

public:
  virtual Types::Coordinate GetParamStep( const size_t idx, const Types::Coordinate mmStep = 1 ) const 
  {
    return Warp->GetParamStep( idx, Vector3D( this->FloatingSize ), mmStep );
  }

  virtual size_t ParamVectorDim() const 
  {
    return Warp->ParamVectorDim();
  }

  virtual size_t VariableParamVectorDim() const 
  {
    return Warp->VariableParamVectorDim();
  }

  virtual void SetWarpXform ( typename W::SmartPtr& warp );

  virtual void GetParamVector ( CoordinateVector& v ) 
  {
    Warp->GetParamVector( v );
  }

protected:
  VoxelMatchingElasticFunctional_WarpTemplate( UniformVolume::SmartPtr& reference, UniformVolume::SmartPtr& floating )
    : VoxelMatchingElasticFunctional( reference, floating ),
      Warp( NULL ), InverseTransformation( NULL )
  {}
  
  virtual ~VoxelMatchingElasticFunctional_WarpTemplate() {}
};

template<class VM> 
class VoxelMatchingElasticFunctional_Template 
  : public VoxelMatchingFunctional_Template<VM>,
    public VoxelMatchingElasticFunctional_WarpTemplate<SplineWarpXform> 
{
public:
  typedef VoxelMatchingElasticFunctional_Template<VM> Self;

  VoxelMatchingElasticFunctional_Template( UniformVolume::SmartPtr& reference, UniformVolume::SmartPtr& floating )
    : VoxelMatchingFunctional_Template<VM>( reference, floating ), 
      VoxelMatchingElasticFunctional_WarpTemplate<SplineWarpXform>( reference, floating ),
      m_ForceOutsideFlag( false ),
      m_ForceOutsideValueRescaled( 0 )
  {
    IncrementalMetric = typename VM::SmartPtr( new VM( *(this->Metric) ) );

    WarpedVolume = NULL;

    DimsX = this->ReferenceGrid->GetDims()[0];
    DimsY = this->ReferenceGrid->GetDims()[1];
    DimsZ = this->ReferenceGrid->GetDims()[2];

    FltDimsX = this->FloatingGrid->GetDims()[0];
    FltDimsY = this->FloatingGrid->GetDims()[1];
  }

  virtual ~VoxelMatchingElasticFunctional_Template() 
  {
    if ( WarpedVolume ) 
      Memory::DeleteArray( WarpedVolume );
  }

  virtual void SetForceOutside
  ( const bool flag = true, const Types::DataItem value = 0 )
  {
    this->m_ForceOutsideFlag = flag;
    this->m_ForceOutsideValueRescaled = this->Metric->DataY.ValueToIndex( value );
  }

  typename Self::ReturnType EvaluateIncremental( const SplineWarpXform* warp, SmartPointer<VM>& localMetric, const DataGrid::RegionType& voi ) 
  {
    Vector3D *pVec;
    int pX, pY, pZ, offset, r;
    int fltIdx[3];
    Types::Coordinate fltFrac[3];

    int endLineIncrement = ( voi.From()[0] + (DimsX - voi.To()[0]) );
    int endPlaneIncrement = DimsX * ( voi.From()[1] + (DimsY - voi.To()[2]) );
    
    const typename VM::Exchange unsetY = this->Metric->DataY.padding();
    *localMetric = *this->Metric;
    r = voi.From()[0] + DimsX * ( voi.From()[1] + DimsY * voi.From()[2] );
    for ( pZ = voi.From()[2]; pZ<voi.To()[2]; ++pZ ) 
      {
      for ( pY = voi.From()[1]; pY<voi.To()[1]; ++pY ) 
        {
        pVec = this->VectorCache;
        warp->GetTransformedGridRow( voi.From()[0]-voi.To()[0], pVec, voi.From()[0], pY, pZ );
        for ( pX = voi.From()[0]; pX<voi.To()[0]; ++pX, ++r, ++pVec ) 
          {
          // Remove this sample from incremental metric according to "ground warp" image.
          const typename VM::Exchange sampleX = this->Metric->GetSampleX( r );
          if ( this->WarpedVolume[r] != unsetY )
            localMetric->Decrement( sampleX, WarpedVolume[r] );
            
          // Tell us whether the current location is still within the floating volume and get the respective voxel.
          *pVec *= this->FloatingInverseDelta;
          if ( this->FloatingGrid->FindVoxelByIndex( *pVec, fltIdx, fltFrac ) ) 
            {
            // Compute data index of the floating voxel in the floating volume.
            offset = fltIdx[0] + FltDimsX * ( fltIdx[1] + FltDimsY * fltIdx[2] );
            
            // Continue metric computation.
            localMetric->Increment( sampleX, this->Metric->GetSampleY( offset, fltFrac ) );
            } 
          else 
            {
            if ( this->m_ForceOutsideFlag )
              {
              localMetric->Increment( sampleX, this->m_ForceOutsideValueRescaled );
              }
            }
          }
        r += endLineIncrement;
        }
      r += endPlaneIncrement;
      }
    
    return localMetric->Get();
  }

  void UpdateWarpFixedParameters();

  virtual typename Self::ReturnType EvaluateWithGradient( CoordinateVector& v, CoordinateVector& g, const Types::Coordinate step = 1 ) 
  {
    const typename Self::ReturnType current = this->EvaluateAt( v );

    if ( this->m_AdaptiveFixParameters && this->WarpNeedsFixUpdate ) 
      {
      this->UpdateWarpFixedParameters();
      }
    
    Types::Coordinate *p = this->Warp->m_Parameters;

    Types::Coordinate pOld;
    double upper, lower;
    const DataGrid::RegionType *voi = this->VolumeOfInfluence;
    for ( size_t dim = 0; dim < this->Dim; ++dim, ++voi ) 
      {
      if ( this->StepScaleVector[dim] <= 0 ) 
        {
        g[dim] = 0;
        } 
      else
        {
        pOld = p[dim];
        
        Types::Coordinate thisStep = step * this->StepScaleVector[dim];
        
        p[dim] += thisStep;
        upper = EvaluateIncremental( this->Warp, IncrementalMetric, *voi );
        
        p[dim] = pOld - thisStep;
        lower = EvaluateIncremental( this->Warp, IncrementalMetric, *voi );
        
        p[dim] = pOld;
        this->WeightedDerivative( lower, upper, *(this->Warp), dim, thisStep );
        
        if ( (upper > current ) || (lower > current) ) 
          {
          g[dim] = upper-lower;
          } 
        else 
          {
          g[dim] = 0;
          }
        }
      }
    
    return current;
  }
  
  virtual typename Self::ReturnType EvaluateAt( CoordinateVector& v )
  {
    this->Warp->SetParamVector( v );
    return this->Evaluate();
  }

  virtual typename Self::ReturnType Evaluate()
  {
    if ( ! WarpedVolume ) 
      WarpedVolume = Memory::AllocateArray<typename VM::Exchange>(  DimsX * DimsY * DimsZ  );

    this->Metric->Reset();
    const typename VM::Exchange unsetY = this->Metric->DataY.padding();

    Vector3D *pVec;
    int pX, pY, pZ;
    int fltIdx[3];
    Types::Coordinate fltFrac[3];

    int r = 0;
    for ( pZ = 0; pZ<DimsZ; ++pZ ) 
      {
      for ( pY = 0; pY<DimsY; ++pY )
        {
        this->Warp->GetTransformedGridRow( DimsX, this->VectorCache, 0, pY, pZ );
        pVec = this->VectorCache;
        for ( pX = 0; pX<DimsX; ++pX, ++r, ++pVec )
          {
          // Tell us whether the current location is still within the floating volume and get the respective voxel.
          *pVec *= this->FloatingInverseDelta;
          if ( this->FloatingGrid->FindVoxelByIndex( *pVec, fltIdx, fltFrac ) ) 
            {
            // Compute data index of the floating voxel in the
            // floating volume.
            const size_t offset = fltIdx[0] + FltDimsX * ( fltIdx[1] + FltDimsY * fltIdx[2] );
            
            // Continue metric computation.
            this->WarpedVolume[r] = this->Metric->GetSampleY( offset, fltFrac );
            this->Metric->Increment( this->Metric->GetSampleX( r ), this->WarpedVolume[r] );
            }
          else 
            {
            if ( this->m_ForceOutsideFlag )
              {
              this->WarpedVolume[r] = this->m_ForceOutsideValueRescaled;
              this->Metric->Increment( this->Metric->GetSampleX( r ), this->WarpedVolume[r] );
              }
            else
              {
              this->WarpedVolume[r] = unsetY;
              }
            }
          }
        }
      }
    
    return this->WeightedTotal( this->Metric->Get(), this->Warp );
  }

protected:
  typename VM::Exchange *WarpedVolume;

  bool m_ForceOutsideFlag;

  typename VM::Exchange m_ForceOutsideValueRescaled;

   SmartPointer<VM> IncrementalMetric;
   
   DataGrid::IndexType::ValueType DimsX, DimsY, DimsZ;
   
   DataGrid::IndexType::ValueType FltDimsX, FltDimsY;
};

VoxelMatchingElasticFunctional* CreateElasticFunctional( const int metric, UniformVolume::SmartPtr& refVolume, UniformVolume::SmartPtr& fltVolume );

} // namespace

#endif // __cmtkVoxelMatchingElasticFunctional_h_included_