text: /*=========================================================================
   *
   *  Copyright Insight Software Consortium
   *
   *  Licensed under the Apache License, Version 2.0 (the "License");
   *  you may not use this file except in compliance with the License.
   *  You may obtain a copy of the License at
   *
   *         http://www.apache.org/licenses/LICENSE-2.0.txt
   *
   *  Unless required by applicable law or agreed to in writing, software
   *  distributed under the License is distributed on an "AS IS" BASIS,
   *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   *  See the License for the specific language governing permissions and
   *  limitations under the License.
   *
   *=========================================================================*/
  #ifndef itkBSplineTransformParametersAdaptor_hxx
  #define itkBSplineTransformParametersAdaptor_hxx
  
  #include "itkBSplineTransformParametersAdaptor.h"
  
  #include "itkBSplineDecompositionImageFilter.h"
  #include "itkBSplineResampleImageFunction.h"
  #include "itkResampleImageFilter.h"
  
  namespace itk
  {
  
  template<typename TTransform>
  BSplineTransformParametersAdaptor<TTransform>
  ::BSplineTransformParametersAdaptor()
  {
    this->m_RequiredFixedParameters.SetSize( SpaceDimension * ( SpaceDimension + 3 ) );
  
    this->m_RequiredTransformDomainOrigin.Fill( 0.0 );
    this->m_RequiredTransformDomainDirection.SetIdentity();
    this->m_RequiredTransformDomainPhysicalDimensions.Fill( 1.0 );
    this->m_RequiredTransformDomainMeshSize.Fill( 1 );
  
    this->UpdateRequiredFixedParameters();
  }
  
  template<typename TTransform>
  BSplineTransformParametersAdaptor<TTransform>
  ::~BSplineTransformParametersAdaptor()
  {
  }
  
  template<typename TTransform>
  void
  BSplineTransformParametersAdaptor<TTransform>
  ::SetRequiredTransformDomainOrigin( const OriginType & origin )
  {
    if( origin != this->m_RequiredTransformDomainOrigin )
      {
      itkDebugMacro( "Setting m_RequiredTransformDomainOrigin to " << origin );
      this->m_RequiredTransformDomainOrigin = origin;
      this->UpdateRequiredFixedParameters();
  
      this->Modified();
      }
  }
  
  template<typename TTransform>
  void
  BSplineTransformParametersAdaptor<TTransform>
  ::SetRequiredTransformDomainPhysicalDimensions( const PhysicalDimensionsType & dimensions )
  {
    if( dimensions != this->m_RequiredTransformDomainPhysicalDimensions )
      {
      itkDebugMacro( "Setting m_RequiredTransformDomainPhysicalDimensions to " << dimensions );
      this->m_RequiredTransformDomainPhysicalDimensions = dimensions;
      this->UpdateRequiredFixedParameters();
  
      this->Modified();
      }
  }
  
  template<typename TTransform>
  void
  BSplineTransformParametersAdaptor<TTransform>
  ::SetRequiredTransformDomainDirection( const DirectionType & direction )
  {
    if( direction != this->m_RequiredTransformDomainDirection )
      {
      itkDebugMacro( "Setting m_RequiredTransformDomainDirection to " << direction );
      this->m_RequiredTransformDomainDirection = direction;
      this->UpdateRequiredFixedParameters();
  
      this->Modified();
      }
  }
  
  template<typename TTransform>
  void
  BSplineTransformParametersAdaptor<TTransform>
  ::SetRequiredTransformDomainMeshSize( const MeshSizeType & meshSize )
  {
    if( meshSize != this->m_RequiredTransformDomainMeshSize )
      {
      itkDebugMacro( "Setting m_RequiredTransformDomainMeshSize to " << meshSize );
      this->m_RequiredTransformDomainMeshSize = meshSize;
      this->UpdateRequiredFixedParameters();
  
      this->Modified();
      }
  }
  
  template<typename TTransform>
  void
  BSplineTransformParametersAdaptor<TTransform>
  ::SetRequiredFixedParameters( const ParametersType fixedParameters )
  {
    Superclass::SetRequiredFixedParameters( fixedParameters );
  
    // Set the direction parameters
    for( SizeValueType di = 0; di < SpaceDimension; di++ )
      {
      for( SizeValueType dj = 0; dj < SpaceDimension; dj++ )
        {
        this->m_RequiredTransformDomainDirection[di][dj] =
          this->m_RequiredFixedParameters[3 * SpaceDimension + ( di * SpaceDimension + dj )];
        }
      }
  
    // set the mesh size parameters
    for( SizeValueType i = 0; i < SpaceDimension; i++ )
      {
      this->m_RequiredTransformDomainMeshSize[i] = static_cast<SizeValueType>( this->m_RequiredFixedParameters[i] ) - TransformType::SplineOrder;
      }
  
    // Set the physical dimensions parameters
    for( SizeValueType i = 0; i < SpaceDimension; i++ )
      {
      ParametersValueType gridSpacing = this->m_RequiredFixedParameters[2 * SpaceDimension + i];
      this->m_RequiredTransformDomainPhysicalDimensions[i] = gridSpacing *
        static_cast<ParametersValueType>( this->m_RequiredTransformDomainMeshSize[i] );
      }
  
    // Set the origin parameters
    OriginType origin;
    for( SizeValueType i = 0; i < SpaceDimension; i++ )
      {
      ParametersValueType gridSpacing = this->m_RequiredFixedParameters[2 * SpaceDimension + i];
      origin[i] = 0.5 * gridSpacing * ( TransformType::SplineOrder - 1 );
      }
    origin = this->m_RequiredTransformDomainDirection * origin;
    for( SizeValueType i = 0; i < SpaceDimension; i++ )
      {
      this->m_RequiredTransformDomainOrigin[i] = origin[i] + this->m_RequiredFixedParameters[SpaceDimension + i];
      }
  }
  
  
  template<typename TTransform>
  void
  BSplineTransformParametersAdaptor<TTransform>
  ::UpdateRequiredFixedParameters()
  {
    // Fixed parameters store the following information:
    //  grid size
    //  grid origin
    //  grid spacing
    //  grid direction
    //  The size of these is equal to the SpaceDimension
  
    // set the grid size parameters
    for( SizeValueType i = 0; i < SpaceDimension; i++ )
      {
      this->m_RequiredFixedParameters[i] = this->m_RequiredTransformDomainMeshSize[i] + TransformType::SplineOrder;
      }
  
    // Set the origin parameters
    OriginType origin;
    for( SizeValueType i = 0; i < SpaceDimension; i++ )
      {
      ParametersValueType gridSpacing = this->m_RequiredTransformDomainPhysicalDimensions[i] /
        static_cast<ParametersValueType>( this->m_RequiredTransformDomainMeshSize[i] );
      origin[i] = -0.5 * gridSpacing * ( TransformType::SplineOrder - 1 );
      }
    origin = this->m_RequiredTransformDomainDirection * origin;
    for( SizeValueType i = 0; i < SpaceDimension; i++ )
      {
      this->m_RequiredFixedParameters[SpaceDimension + i] = static_cast<ParametersValueType>(
        origin[i] + this->m_RequiredTransformDomainOrigin[i] );
      }
  
    // Set the spacing parameters
    for( SizeValueType i = 0; i < SpaceDimension; i++ )
      {
      ParametersValueType gridSpacing = this->m_RequiredTransformDomainPhysicalDimensions[i] /
        static_cast<ParametersValueType>( this->m_RequiredTransformDomainMeshSize[i] );
      this->m_RequiredFixedParameters[2 * SpaceDimension + i] =
        static_cast<ParametersValueType>( gridSpacing );
      }
  
    // Set the direction parameters
    for( SizeValueType di = 0; di < SpaceDimension; di++ )
      {
      for( SizeValueType dj = 0; dj < SpaceDimension; dj++ )
        {
        this->m_RequiredFixedParameters[3 * SpaceDimension + ( di * SpaceDimension + dj )] =
          static_cast<ParametersValueType>( this->m_RequiredTransformDomainDirection[di][dj] );
        }
      }
  }
  
  template<typename TTransform>
  void
  BSplineTransformParametersAdaptor<TTransform>
  ::AdaptTransformParameters()
  {
    if( !this->m_Transform )
      {
      itkExceptionMacro( "Transform has not been set." );
      return;
      }
  
    if( this->m_RequiredFixedParameters == this->m_Transform->GetFixedParameters() )
      {
      return;
      }
  
    SizeType newGridSize;
    OriginType newGridOrigin;
    SpacingType newGridSpacing;
    DirectionType newGridDirection;
    for( SizeValueType i = 0; i < SpaceDimension; i++ )
      {
      newGridSize[i] = static_cast<SizeValueType>( this->m_RequiredFixedParameters[i] );
      newGridOrigin[i] = this->m_RequiredFixedParameters[SpaceDimension + i];
      newGridSpacing[i] = this->m_RequiredFixedParameters[2 * SpaceDimension + i];
      for( SizeValueType j = 0; j < SpaceDimension; j++ )
        {
        newGridDirection[i][j] = this->m_RequiredFixedParameters[3 * SpaceDimension + ( i * SpaceDimension + j )];
        }
      }
  
    const RegionType & coefficientImageRegion =
      this->m_Transform->GetCoefficientImages()[0]->GetLargestPossibleRegion();
    IndexType newGridIndex = coefficientImageRegion.GetIndex();
  
    //  Resample the coefficient images
    typedef BSplineResampleImageFunction<ImageType, ParametersValueType> CoefficientUpsampleFunctionType;
    typedef ResampleImageFilter<ImageType, ImageType, ParametersValueType> UpsampleFilterType;
    typedef BSplineDecompositionImageFilter<ImageType, ImageType> DecompositionFilterType;
  
    CoefficientImageArray newCoefficientImages;
  
    // Loop over dimension: each direction is upsampled separately.
    for( SizeValueType j = 0; j < SpaceDimension; j++ )
      {
      /* Set the coeficient image as the input of the upsampler filter.
       * The upsampler samples the deformation field at the locations
       * of the new control points, given the current coefficients
       * (note: it does not just interpolate the coefficient image,
       * which would be wrong). The B-spline coefficients that
       * describe the resulting image are computed by the
       * decomposition filter.
       *
       * This code is copied from the itk-example
       * DeformableRegistration6.cxx .
       */
      typename UpsampleFilterType::Pointer upsampler = UpsampleFilterType::New();
      typename CoefficientUpsampleFunctionType::Pointer coeffUpsampleFunction =
        CoefficientUpsampleFunctionType::New();
  
      upsampler->SetInterpolator( coeffUpsampleFunction );
      upsampler->SetSize( newGridSize );
      upsampler->SetOutputStartIndex( newGridIndex );
      upsampler->SetOutputSpacing( newGridSpacing );
      upsampler->SetOutputOrigin( newGridOrigin );
      upsampler->SetOutputDirection( newGridDirection );
      upsampler->SetInput( this->m_Transform->GetCoefficientImages()[j] );
  
      typename DecompositionFilterType::Pointer decompositionFilter = DecompositionFilterType::New();
      decompositionFilter->SetSplineOrder( TransformType::SplineOrder );
      decompositionFilter->SetInput( upsampler->GetOutput() );
  
      newCoefficientImages[j] = decompositionFilter->GetOutput();
      newCoefficientImages[j]->Update();
      newCoefficientImages[j]->DisconnectPipeline();
      }
  
    this->m_Transform->SetCoefficientImages( newCoefficientImages );
  }
  
  template<typename TTransform>
  void
  BSplineTransformParametersAdaptor<TTransform>
  ::PrintSelf( std::ostream & os, Indent indent ) const
  {
    Superclass::PrintSelf( os, indent );
  
    os << "Required transform domain origin: " << this->m_RequiredTransformDomainOrigin << std::endl;
    os << "Required transform domain direction: " << this->m_RequiredTransformDomainDirection << std::endl;
    os << "Required transform domain physical dimensions: " << this->m_RequiredTransformDomainPhysicalDimensions << std::endl;
    os << "Required transform domain mesh size: " << this->m_RequiredTransformDomainMeshSize << std::endl;
  }
  
  }  // namespace itk
  
  #endif
  