cmtkAffineXform.h

Go to the documentation of this file.

/*
//
//  Copyright 1997-2009 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 __cmtkAffineXform_h_included_
#define __cmtkAffineXform_h_included_

#include <cmtkconfig.h>

#include <Base/cmtkTypes.h>
#include <Base/cmtkVector.h>
#include <Base/cmtkXform.h>
#include <Base/cmtkMatrix4x4.h>

#include <System/cmtkSmartPtr.h>

#include <Base/cmtkUnits.h>

namespace
cmtk
{


class AffineXform : 
  public Xform 
{
public:
  typedef AffineXform Self;

  typedef Xform Superclass;

  typedef SmartPointer<Self> SmartPtr;

  typedef SmartConstPointer<Self> SmartConstPtr;

  typedef Matrix4x4<Types::Coordinate> MatrixType;

  void MakeIdentityXform ();

  Self::MatrixType Matrix;

  static const size_t TotalNumberOfParameters = 15;

  Self::SpaceVectorType RotateScaleShear ( const Self::SpaceVectorType& v ) const;

  AffineXform () :
    m_LogScaleFactors( false )
  {
    this->AllocateParameterVector( TotalNumberOfParameters );
    this->NumberDOFs = this->DefaultNumberOfDOFs();
    this->MakeIdentityXform();
  }

  AffineXform ( const CoordinateVector& v , 
                const bool logScaleFactors = false  ) :
    m_LogScaleFactors( logScaleFactors )
  {
    this->AllocateParameterVector( TotalNumberOfParameters );
    this->NumberDOFs = this->DefaultNumberOfDOFs();
    this->SetParamVector( v );
  }

  AffineXform ( const Types::Coordinate v[15] , 
                const bool logScaleFactors = false  ) :
    m_LogScaleFactors( logScaleFactors )
  {
    this->AllocateParameterVector( TotalNumberOfParameters );
    this->NumberDOFs = this->DefaultNumberOfDOFs();
    memcpy( this->m_Parameters, v, 15 * sizeof(Types::Coordinate) );
    this->ComposeMatrix();
    this->CanonicalRotationRange();
  }

  AffineXform ( const Types::Coordinate matrix[4][4], const Types::Coordinate* center = NULL );

  AffineXform ( const MatrixType& matrix, const Types::Coordinate* center = NULL );

  AffineXform ( const Types::Coordinate matrix[4][4], const Types::Coordinate xlate[3], const Types::Coordinate center[3] );
  
  AffineXform ( const AffineXform& other );
  
  virtual ~AffineXform() 
  { 
    InverseXform = Self::SmartPtr::Null; 
  }
  
  Self::SmartPtr Clone () const 
  {
    return Self::SmartPtr( this->CloneVirtual() );
  }

  virtual Self* MakeInverse () const;

  virtual Self::SmartPtr GetDifference( const AffineXform& other ) const;

  const Self::SmartPtr GetInverse() const;

  virtual Types::Coordinate GetGlobalScaling() const 
  { 
    if ( this->m_LogScaleFactors )
      {
      return exp( this->m_Parameters[6] + this->m_Parameters[7] + this->m_Parameters[8] );
      }
    else
      {
      return this->m_Parameters[6] * this->m_Parameters[7] * this->m_Parameters[8];
      }
  }

  virtual Types::Coordinate GetJacobianDeterminant ( const Self::SpaceVectorType& ) const
  { 
    return this->GetGlobalScaling(); 
  }

  void Concat( const AffineXform& other );

  void Insert( const AffineXform& other );

  void RotateWXYZ( const Units::Radians angle, const Self::SpaceVectorType& direction, const Types::Coordinate* origin = NULL, Self::MatrixType *const accumulate = NULL );

  void ChangeCoordinateSystem
  ( const Self::SpaceVectorType& newX, const Self::SpaceVectorType& newY )
  {
    this->Matrix.ChangeCoordinateSystem( newX, newY );
    this->DecomposeMatrix();
  }

  virtual Self::SpaceVectorType Apply ( const Self::SpaceVectorType& vec ) const 
  {
    Self::SpaceVectorType Result( vec );
    this->ApplyInPlace( Result );
    return Result;
  }
  
  void ApplyInPlace( Self::SpaceVectorType& vec ) const 
  {
    vec *= this->Matrix;
  }

  virtual bool ApplyInverse ( const Self::SpaceVectorType& v, Self::SpaceVectorType& u, const Types::Coordinate = 0.01  ) const
  {
    u = this->GetInverse()->Apply( v );
    return true;
  }

  virtual bool ApplyInverseInPlace( Self::SpaceVectorType& v, const Types::Coordinate = 0.01  ) const
  {
    this->GetInverse()->ApplyInPlace( v );
    return true;
  }


  const Types::Coordinate* RetXlate () const  { return this->m_Parameters; }
  const Types::Coordinate* RetAngles () const { return this->m_Parameters+3; }
  const Types::Coordinate* RetScales () const { return this->m_Parameters+6; }
  const Types::Coordinate* RetShears () const { return this->m_Parameters+9; }
  const Types::Coordinate* RetCenter () const { return this->m_Parameters+12; }


  Types::Coordinate* RetXlate ()  { return this->m_Parameters; }
  Types::Coordinate* RetAngles () { return this->m_Parameters+3; }
  Types::Coordinate* RetScales () { return this->m_Parameters+6; }
  Types::Coordinate* RetShears () { return this->m_Parameters+9; }
  Types::Coordinate* RetCenter () { return this->m_Parameters+12; }


  void SetXlate ( const Types::Coordinate* xlate ) 
  { 
    if ( this->RetXlate() != xlate )
      memcpy( this->RetXlate(), xlate, 3 * sizeof(Types::Coordinate) );
    this->ComposeMatrix();
  }

  void SetXlate ( const Types::Coordinate dx, const Types::Coordinate dy, const Types::Coordinate dz ) 
  { 
    this->m_Parameters[0] = dx; this->m_Parameters[1] = dy; this->m_Parameters[2] = dz;
    this->ComposeMatrix(); 
  }

  void Translate( const Types::Coordinate dx, const Types::Coordinate dy, const Types::Coordinate dz ) 
  {
    this->m_Parameters[0] += dx; this->m_Parameters[1] += dy; this->m_Parameters[2] += dz; 
    this->ComposeMatrix(); 
  }
  
  void Translate( const Self::SpaceVectorType& delta ) 
  { 
    for ( int dim = 0; dim < 3; ++dim ) 
      this->m_Parameters[dim] += delta[dim];
    this->ComposeMatrix();
  }

  void SetAngles ( const Types::Coordinate* angles ) 
  {
    if ( angles != this->RetAngles() )
      memcpy( this->RetAngles(), angles, 3 * sizeof(Types::Coordinate) );
    this->ComposeMatrix();
  }

  template<class T> void GetScales( T (&scales)[3] ) const
  { 
    if ( this->m_LogScaleFactors )
      {
      for ( size_t i = 0; i < 3; ++i )
        {
        scales[i] = exp( this->m_Parameters[6+i] );
        }
      }
    else
      {
      for ( size_t i = 0; i < 3; ++i )
        {
        scales[i] = this->m_Parameters[6+i];
        }
      }
  }

  void SetScales ( const Types::Coordinate* scales ) 
  { 
    if ( this->RetScales() != scales )
      memcpy( this->RetScales(), scales, 3 * sizeof(Types::Coordinate) );
    this->ComposeMatrix();
  }
  
  void SetScales ( const Types::Coordinate sx, const Types::Coordinate sy, const Types::Coordinate sz )
  { 
    this->m_Parameters[6] = sx; 
    this->m_Parameters[7] = sy; 
    this->m_Parameters[8] = sz; }

  void SetShears ( const Types::Coordinate* shears ) 
  { 
    if ( this->RetShears() != shears )
      memcpy( this->RetShears(), shears, 3 * sizeof(Types::Coordinate) );
    this->ComposeMatrix();
  }

  void SetCenter ( const Types::Coordinate* center ) 
  {
    if ( this->RetCenter() != center )
      memcpy( RetCenter(), center, 3 * sizeof(Types::Coordinate) );
    this->ComposeMatrix();
  }
  

  const Types::Coordinate* RetMatrix () const { return &Matrix[0][0]; }
    
  Types::Coordinate* RetMatrix () { return &Matrix[0][0]; }

  void SetMatrix( const float matrix[4][4] );
  void SetMatrix( const double matrix[4][4] );
  void SetMatrix( const MatrixType& matrix );

  template<class T> void GetMatrix( T (&matrix)[4][4] ) const;

  void ChangeCenter ( const Self::SpaceVectorType& center );

  virtual size_t ParamVectorDim () const { return 15; }

  virtual size_t VariableParamVectorDim () const { return std::min( 12, NumberDOFs ); }

  virtual void SetNumberDOFs ( const int numberDOFs = 12 );

  virtual unsigned int GetNumberDOFs () const { return NumberDOFs; }

  bool GetUseLogScaleFactors() const
  {
    return this->m_LogScaleFactors;
  }

  void SetUseLogScaleFactors( const bool logScaleFactors );

  virtual void SetParamVector ( CoordinateVector& v );

  virtual void SetParamVector ( const CoordinateVector& v );

  virtual void SetParameter ( const size_t idx, const Types::Coordinate p );

  virtual Types::Coordinate GetParamStep( const size_t idx, const Self::SpaceVectorType& volSize, const Types::Coordinate step_mm = 1 ) const;
  
  AffineXform& operator=( const AffineXform& other );

  virtual void Print() const;

protected:
  virtual Self* CloneVirtual() const
  {
    return new AffineXform( *this );
  }
  
  void ComposeMatrix();

  bool DecomposeMatrix();

  int NumberDOFs;

private:
  bool m_LogScaleFactors;

  static int DefaultNumberOfDOFs() { return 12; }

  mutable Self::SmartPtr InverseXform;

  void UpdateInverse() const;

  void SetMatrixDirect ( const Types::Coordinate* matrix ) 
  {
    Matrix.Set( matrix );
    this->DecomposeMatrix();
  }

  void CanonicalRotationRange();
};


} // namespace cmtk

#endif // #define __cmtkAffineXform_h_included_