cmtkSplineWarpXform.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 __cmtkSplineWarpXform_h_included_
#define __cmtkSplineWarpXform_h_included_

#include <cmtkconfig.h>

#include <Base/cmtkWarpXform.h>
#include <Base/cmtkVector.h>
#include <Base/cmtkAffineXform.h>
#include <Base/cmtkCubicSpline.h>

#include <cassert>
#include <vector>
#include <algorithm>
  
#include <System/cmtkSmartPtr.h>
#include <System/cmtkThreads.h>

namespace
cmtk
{


class SplineWarpXform :
  public WarpXform 
{
public:
  typedef SplineWarpXform Self;

  typedef WarpXform Superclass;

  typedef SmartPointer<Self> SmartPtr;

  typedef SmartConstPointer<Self> SmartConstPtr;

  SplineWarpXform();

  SplineWarpXform( const FixedVector<3,Types::Coordinate>& domain, const Types::Coordinate delta, const AffineXform *initialXform = NULL, const bool exactDelta = false );

  void Init( const Self::SpaceVectorType& domain, const Types::Coordinate delta, const AffineXform *initialXform = NULL, const bool exactDelta = false );
  
  SplineWarpXform( const FixedVector<3,Types::Coordinate>& domain, const Self::IndexType& dims, CoordinateVector::SmartPtr& parameters, const AffineXform *initialXform = NULL );

  Self::SmartPtr Clone () const 
  {
    return Self::SmartPtr( this->CloneVirtual() );
  }

  void InitControlPoints( const AffineXform* affineXform = NULL );

  virtual void Update( const bool exactDelta = false );

  virtual SplineWarpXform* MakeInverse() const 
  {
    return NULL;
  }
  
  virtual void Refine();

  virtual Types::Coordinate GetGridEnergy() const;

  virtual Types::Coordinate GetGridEnergy( const Types::Coordinate *cp ) const;

  virtual Types::Coordinate GetGridEnergy( const Self::SpaceVectorType& v ) const;

  virtual void GetGridEnergyDerivative( double& lower, double& upper, const int param, const Types::Coordinate step ) const;

  virtual Types::Coordinate GetJacobianDeterminant ( const Self::SpaceVectorType& v ) const;

  virtual Types::Coordinate GetJacobianDeterminant ( const int x, const int y, const int z ) const;

  virtual void GetJacobianDeterminantRow( double *const values, const int x, const int y, const int z, const size_t numberOfPoints = 1 ) const;

  Types::Coordinate JacobianDeterminant ( const Types::Coordinate *cp ) const;

  virtual Types::Coordinate GetJacobianConstraint() const;

  virtual Types::Coordinate GetJacobianFoldingConstraint() const;

  virtual void RelaxToUnfold();

  virtual Types::Coordinate GetRigidityConstraint() const;

  virtual Types::Coordinate GetRigidityConstraint( const DataGrid* weightMap ) const;

  virtual Types::Coordinate GetJacobianConstraintSparse() const;

  virtual Types::Coordinate GetRigidityConstraintSparse() const;

  virtual void GetJacobianConstraintDerivative( double& lower, double& upper, const int param, const UniformVolume::RegionType&, const Types::Coordinate step ) const;

  virtual void GetJacobianFoldingConstraintDerivative( double& lower, double& upper, const int param, const UniformVolume::RegionType&, const Types::Coordinate step ) const;

  virtual void GetJacobianConstraintDerivative( double& lower, double& upper, const int param, const Types::Coordinate step ) const;
  
  virtual void GetRigidityConstraintDerivative( double& lower, double& upper, const int param, const UniformVolume::RegionType&, const Types::Coordinate step ) const;
  
  virtual void GetRigidityConstraintDerivative( double& lower, double& upper, const int param, const UniformVolume::RegionType&, const Types::Coordinate step,
                                                const DataGrid* weightMap ) const;

  virtual void GetRigidityConstraintDerivative( double& lower, double& upper, const int param, const Types::Coordinate step ) const;

  virtual Types::Coordinate GetInverseConsistencyError( const WarpXform* inverse, const UniformVolume* volume, const UniformVolume::RegionType* voi = NULL ) const;
  
  virtual bool ApplyInverse ( const Self::SpaceVectorType& v, Self::SpaceVectorType& u, const Types::Coordinate accuracy = 0.01  ) const;

  virtual bool ApplyInverseInPlace( Self::SpaceVectorType& v, const Types::Coordinate accuracy = 0.01  ) const;

  virtual bool ApplyInverseInPlaceWithInitial( Self::SpaceVectorType& v, const Self::SpaceVectorType& initial, const Types::Coordinate accuracy = 0.01 ) const;

  virtual void ApplyInPlace( Self::SpaceVectorType& v ) const 
  {
    Types::Coordinate r[3], f[3];
    int grid[3];
    
    // Do some precomputations.
    { 
    for ( int dim = 0; dim<3; ++dim ) 
      {
      // This is the (real-valued) index of the control point grid cell the
      // given location is in.
      r[dim] = this->InverseSpacing[dim] * v[dim];
      // This is the actual cell index.
      grid[dim] = std::min<int>( static_cast<int>( r[dim] ), this->m_Dims[dim]-4 );
      // And here's the relative position within the cell.
      f[dim] = r[dim] - grid[dim];
      }
    }

    // Create a pointer to the front-lower-left corner of the c.p.g. cell.
    const Types::Coordinate* coeff = this->m_Parameters + 3 * ( grid[0] + this->m_Dims[0] * (grid[1] + this->m_Dims[1] * grid[2]) );

    for ( int dim = 0; dim<3; ++dim ) 
      {
      Types::Coordinate mm = 0;
      const Types::Coordinate *coeff_mm = coeff;
      
      // Loop over 4 c.p.g. planes in z-direction.
      for ( int m = 0; m < 4; ++m ) 
        {
        Types::Coordinate ll = 0;
        const Types::Coordinate *coeff_ll = coeff_mm;
        
        // Loop over 4 c.p.g. planes in y-direction.
        for ( int l = 0; l < 4; ++l ) 
          {
          Types::Coordinate kk = 0;
          const Types::Coordinate *coeff_kk = coeff_ll;
          
          // Loop over 4 c.p.g. planes in x-direction.
          for ( int k = 0; k < 4; ++k, coeff_kk+=3 ) 
            {
            kk += CubicSpline::ApproxSpline( k, f[0] ) * (*coeff_kk);
            }
          ll += CubicSpline::ApproxSpline( l, f[1] ) * kk;
          coeff_ll += nextJ;
          }     
        mm += CubicSpline::ApproxSpline( m, f[2] ) * ll;
        coeff_mm += nextK;
        }
      v[dim] = mm;
      ++coeff;
      }
  }

  virtual void GetVolumeOfInfluence( const size_t idx, const Self::SpaceVectorType&, const Self::SpaceVectorType&, Self::SpaceVectorType&, Self::SpaceVectorType&, const int = -1 ) const;
  
  void RegisterVolume( const UniformVolume* volume )
  {
    this->RegisterVolumePoints( volume->m_Dims, volume->m_Delta, volume->m_Offset );
  }

  void UnRegisterVolume();
  
  void GetTransformedGrid( Self::SpaceVectorType& v, const int idxX, const int idxY, const int idxZ ) const;
  
  void GetTransformedGridRow( const int numPoints, Self::SpaceVectorType *const v, const int idxX, const int idxY, const int idxZ ) const;
  
  virtual Types::Coordinate GetParamStep( const size_t idx, const Self::SpaceVectorType& volSize, const Types::Coordinate mmStep = 1 ) const 
  {
    return 4 * WarpXform::GetParamStep( idx, volSize, mmStep );
  }
  
  virtual Self::SpaceVectorType& GetDeformedControlPointPosition( Self::SpaceVectorType&, const int, const int, const int ) const;
  
  Types::Coordinate* GetPureDeformation( const bool includeScale = false ) const;

  virtual CoordinateMatrix3x3 GetJacobian( const Self::SpaceVectorType& v ) const;

  virtual void GetJacobian( const Self::SpaceVectorType& v, CoordinateMatrix3x3& J ) const;

  virtual void GetJacobianAtControlPoint( const Types::Coordinate* cp, CoordinateMatrix3x3& J ) const;

  virtual void GetJacobianRow( CoordinateMatrix3x3 *const array, const int x, const int y, const int z, const size_t numberOfPoints ) const;
  
private:
  void RegisterVolumePoints ( const DataGrid::IndexType&, const Self::SpaceVectorType& );

  void RegisterVolumePoints( const DataGrid::IndexType&, const Self::SpaceVectorType&, const Self::SpaceVectorType& );

  void RegisterVolumeAxis ( const DataGrid::IndexType::ValueType, const Types::Coordinate delta, const Types::Coordinate origin, const int, const Types::Coordinate, 
                            std::vector<int>& g, std::vector<Types::Coordinate>& spline, std::vector<Types::Coordinate>& dspline );

  Types::Coordinate GetRigidityConstraint( const CoordinateMatrix3x3& J ) const;

protected:
  virtual SplineWarpXform* CloneVirtual () const;

  DataGrid::IndexType VolumeDims;


  std::vector<int> gX;
  std::vector<int> gY;
  std::vector<int> gZ;


  std::vector<Types::Coordinate> splineX;
  std::vector<Types::Coordinate> splineY;
  std::vector<Types::Coordinate> splineZ;


  std::vector<Types::Coordinate> dsplineX;
  std::vector<Types::Coordinate> dsplineY;
  std::vector<Types::Coordinate> dsplineZ;

  int GridPointOffset[48];

  void Init ();

  typedef struct 
  {
    const SplineWarpXform *thisObject;
    int ThisThreadIndex;
    int NumberOfThreads;
    Types::Coordinate Constraint;
  } JacobianConstraintThreadInfo;
  
  static void GetJacobianConstraintThread( void *const args, const size_t taskIdx, const size_t taskCnt, const size_t, const size_t );

  static void GetJacobianFoldingConstraintThread( void *const args, const size_t taskIdx, const size_t taskCnt, const size_t, const size_t );

  void FindClosestControlPoint( const Self::SpaceVectorType& v, Self::SpaceVectorType& cp ) const;

  friend class SplineWarpXformUniformVolume;
};

} // namespace

#endif // #ifndef __cmtkSplineWarpXform_h_included_