itkRecursiveBSplineTransformImplementation.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright UMC Utrecht and contributors
 *
 *  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 __itkRecursiveBSplineTransformImplementation_h
#define __itkRecursiveBSplineTransformImplementation_h
 
#include "itkRecursiveBSplineInterpolationWeightFunction.h"
 
namespace itk
{
 
template< unsigned int OutputDimension, unsigned int SpaceDimension, unsigned int SplineOrder, class TScalar >
class RecursiveBSplineTransformImplementation
{
public:
 
  typedef TScalar ScalarType;
  typedef double  InternalFloatType;
 
  itkStaticConstMacro( HelperConstVariable, unsigned int,
    ( SpaceDimension - 1 ) * ( SplineOrder + 1 ) );
 
  typedef itk::RecursiveBSplineInterpolationWeightFunction<
    ScalarType, OutputDimension, SplineOrder > RecursiveBSplineWeightFunctionType;
  itkStaticConstMacro( BSplineNumberOfIndices, unsigned int,
    RecursiveBSplineWeightFunctionType::NumberOfIndices );
 
  typedef ScalarType *  OutputPointType;
  typedef ScalarType ** CoefficientPointerVectorType;
 
  static inline void TransformPoint(
    OutputPointType opp, const CoefficientPointerVectorType mu,
    const OffsetValueType * gridOffsetTable,
    const double * weights1D )
  {
    ScalarType * tmp_mu[ OutputDimension ];
    for( unsigned int j = 0; j < OutputDimension; ++j )
    {
      tmp_mu[ j ] = mu[ j ];
    }
 
    ScalarType tmp_opp[ OutputDimension ];
    for( unsigned int j = 0; j < OutputDimension; ++j )
    {
      opp[ j ] = 0.0;
    }
 
    OffsetValueType bot = gridOffsetTable[ SpaceDimension - 1 ];
    for( unsigned int k = 0; k <= SplineOrder; ++k )
    {
      RecursiveBSplineTransformImplementation< OutputDimension, SpaceDimension - 1, SplineOrder, TScalar >
        ::TransformPoint( tmp_opp, tmp_mu, gridOffsetTable, weights1D );
 
      for( unsigned int j = 0; j < OutputDimension; ++j )
      {
        opp[ j ] += tmp_opp[ j ] * weights1D[ k + HelperConstVariable ];
 
        // move to the next mu
        tmp_mu[ j ] += bot;
      }
    }
  } // end TransformPoint()
 
 
  static inline void GetJacobian(
    ScalarType * & jacobians, const double * weights1D, double value )
  {
    for( unsigned int k = 0; k <= SplineOrder; ++k )
    {
      RecursiveBSplineTransformImplementation< OutputDimension, SpaceDimension - 1, SplineOrder, TScalar >
        ::GetJacobian( jacobians, weights1D, value * weights1D[ k + HelperConstVariable ] );
    }
  } // end GetJacobian()
 
 
  static inline void EvaluateJacobianWithImageGradientProduct(
    ScalarType * & imageJacobian, const InternalFloatType * movingImageGradient,
    const double * weights1D, double value )
  {
    for( unsigned int k = 0; k <= SplineOrder; ++k )
    {
      RecursiveBSplineTransformImplementation< OutputDimension, SpaceDimension - 1, SplineOrder, TScalar >
        ::EvaluateJacobianWithImageGradientProduct( imageJacobian, movingImageGradient, weights1D,
        value * weights1D[ k + HelperConstVariable ] );
    }
  } // end EvaluateJacobianWithImageGradientProduct()
 
 
  static inline void ComputeNonZeroJacobianIndices(
    unsigned long * & nzji,
    const unsigned long parametersPerDim,
    unsigned long currentIndex,
    const OffsetValueType * gridOffsetTable )
  {
    const OffsetValueType bot = gridOffsetTable[ SpaceDimension - 1 ];
    for( unsigned int k = 0; k <= SplineOrder; ++k )
    {
      RecursiveBSplineTransformImplementation< OutputDimension, SpaceDimension - 1, SplineOrder, TScalar >
        ::ComputeNonZeroJacobianIndices( nzji, parametersPerDim, currentIndex, gridOffsetTable );
 
      currentIndex += bot;
    }
  } // end ComputeNonZeroJacobianIndices()
 
 
  static inline void GetSpatialJacobian(
    InternalFloatType * sj,
    const CoefficientPointerVectorType mu,
    const OffsetValueType * gridOffsetTable,
    const double * weights1D,                    // normal B-spline weights
    const double * derivativeWeights1D )         // 1st derivative of B-spline
  {
    ScalarType * tmp_mu[ OutputDimension ];
    for( unsigned int j = 0; j < OutputDimension; ++j )
    {
      tmp_mu[ j ] = mu[ j ];
    }
 
    InternalFloatType tmp_sj[ OutputDimension * SpaceDimension ];
    for( unsigned int n = 0; n < OutputDimension * ( SpaceDimension + 1 ); ++n )
    {
      sj[ n ] = 0.0;
    }
 
    OffsetValueType bot = gridOffsetTable[ SpaceDimension - 1 ];
    for( unsigned int k = 0; k <= SplineOrder; ++k )
    {
      RecursiveBSplineTransformImplementation< OutputDimension, SpaceDimension - 1, SplineOrder, TScalar >
        ::GetSpatialJacobian( tmp_sj, tmp_mu, gridOffsetTable, weights1D, derivativeWeights1D );
 
      for( unsigned int n = 0; n < OutputDimension * SpaceDimension; ++n )
      {
        sj[ n ] += tmp_sj[ n ] * weights1D[ k + HelperConstVariable ];
      }
 
      for( unsigned int j = 0; j < OutputDimension; ++j )
      {
        sj[ OutputDimension * SpaceDimension + j ]
          += tmp_sj[ j ] * derivativeWeights1D[ k + HelperConstVariable ];
 
        // move to the next mu
        tmp_mu[ j ] += bot;
      }
    }
  } // end GetSpatialJacobian()
 
 
  static inline void GetSpatialHessian(
    InternalFloatType * sh,
    const CoefficientPointerVectorType mu,
    const OffsetValueType * gridOffsetTable,
    const double * weights1D,                   // normal B-spline weights
    const double * derivativeWeights1D,         // 1st derivative of B-spline
    const double * hessianWeights1D )           // 2nd derivative of B-spline
  {
    const unsigned int helperDim1 = OutputDimension * SpaceDimension * ( SpaceDimension + 1 ) / 2;
    const unsigned int helperDim2 = OutputDimension * ( SpaceDimension + 1 ) * ( SpaceDimension + 2 ) / 2;
 
    ScalarType * tmp_mu[ OutputDimension ];
    for( unsigned int j = 0; j < OutputDimension; ++j )
    {
      tmp_mu[ j ] = mu[ j ];
    }
 
    InternalFloatType tmp_sh[ helperDim1 ];
    for( unsigned int n = 0; n < helperDim2; ++n )
    {
      sh[ n ] = 0.0;
    }
 
    OffsetValueType bot = gridOffsetTable[ SpaceDimension - 1 ];
    for( unsigned int k = 0; k <= SplineOrder; ++k )
    {
      RecursiveBSplineTransformImplementation< OutputDimension, SpaceDimension - 1, SplineOrder, TScalar >
        ::GetSpatialHessian( tmp_sh, tmp_mu, gridOffsetTable, weights1D, derivativeWeights1D, hessianWeights1D );
 
      for( unsigned int n = 0; n < helperDim1; ++n )
      {
        sh[ n ] += tmp_sh[ n ] * weights1D[ k + HelperConstVariable ];
      }
 
      for( unsigned int n = 0; n < SpaceDimension; ++n )
      {
        for( unsigned int j = 0; j < OutputDimension; ++j )
        {
          sh[ OutputDimension * n + helperDim1 + j ]
            += tmp_sh[ OutputDimension * n * ( n + 1 ) / 2 + j ] * derivativeWeights1D[ k + HelperConstVariable ];
        }
      }
 
      for( unsigned int j = 0; j < OutputDimension; ++j )
      {
        sh[ helperDim2 - OutputDimension + j ]
          += tmp_sh[ j ] * hessianWeights1D[ k + HelperConstVariable ];
 
        // move to the next mu
        tmp_mu[ j ] += bot;
      }
    }
  } // end GetSpatialHessian()
 
 
  static inline void GetJacobianOfSpatialJacobian(
    InternalFloatType * & jsj_out,
    const double * weights1D,                   // normal B-spline weights
    const double * derivativeWeights1D,         // 1st derivative of B-spline
    const double * directionCosines,
    InternalFloatType * jsj )
  {
    const unsigned int helperDim = OutputDimension - SpaceDimension + 1;
 
    InternalFloatType tmp_jsj[ helperDim + 1 ];
 
    for( unsigned int k = 0; k <= SplineOrder; ++k )
    {
      const double w  = weights1D[ k + HelperConstVariable ];
      const double dw = derivativeWeights1D[ k + HelperConstVariable ];
 
      for( unsigned int n = 0; n < helperDim; ++n )
      {
        tmp_jsj[ n ] = jsj[ n ] * w;
      }
 
      tmp_jsj[ helperDim ] = jsj[ 0 ] * dw;
 
      RecursiveBSplineTransformImplementation< OutputDimension, SpaceDimension - 1, SplineOrder, TScalar >
        ::GetJacobianOfSpatialJacobian( jsj_out, weights1D, derivativeWeights1D, directionCosines, tmp_jsj );
    }
  } // end GetJacobianOfSpatialJacobian()
 
 
  static inline void GetJacobianOfSpatialHessian(
    InternalFloatType * & jsh_out,
    const double * weights1D,                   // normal B-spline weights
    const double * derivativeWeights1D,         // 1st derivative of B-spline
    const double * hessianWeights1D,            // 2nd derivative of B-spline
    const double * directionCosines,
    InternalFloatType * jsh )
  {
    const unsigned int helperDim   = OutputDimension - SpaceDimension;
    const unsigned int helperDimW  = ( helperDim + 1 ) * ( helperDim + 2 ) / 2;
    const unsigned int helperDimDW = helperDim + 1;
 
    InternalFloatType tmp_jsh[ helperDimW + helperDimDW + 1 ];
 
    for( unsigned int k = 0; k <= SplineOrder; ++k )
    {
      const double w  = weights1D[ k + HelperConstVariable ];
      const double dw = derivativeWeights1D[ k + HelperConstVariable ];
      const double hw = hessianWeights1D[ k + HelperConstVariable ];
 
      for( unsigned int n = 0; n < helperDimW; ++n )
      {
        tmp_jsh[ n ] = jsh[ n ] * w;
      }
 
      for( unsigned int n = 0; n < helperDimDW; ++n )
      {
        unsigned int nn = n * ( n + 1 ) / 2;
        tmp_jsh[ n + helperDimW ] = jsh[ nn ] * dw;
      }
 
      tmp_jsh[ helperDimW + helperDimDW ] = jsh[ 0 ] * hw;
 
      RecursiveBSplineTransformImplementation< OutputDimension, SpaceDimension - 1, SplineOrder, TScalar >
        ::GetJacobianOfSpatialHessian( jsh_out, weights1D, derivativeWeights1D, hessianWeights1D, directionCosines, tmp_jsh );
    }
  } // end GetJacobianOfSpatialHessian()
 
 
};
 
 
template< unsigned int OutputDimension, unsigned int SplineOrder, class TScalar >
class RecursiveBSplineTransformImplementation< OutputDimension, 0, SplineOrder, TScalar >
{
public:
 
  typedef TScalar ScalarType;
  typedef double  InternalFloatType;
 
  typedef itk::RecursiveBSplineInterpolationWeightFunction<
    TScalar, OutputDimension, SplineOrder > RecursiveBSplineWeightFunctionType;
  itkStaticConstMacro( BSplineNumberOfIndices, unsigned int,
    RecursiveBSplineWeightFunctionType::NumberOfIndices );
 
  typedef ScalarType *  OutputPointType;
  typedef ScalarType ** CoefficientPointerVectorType;
 
  static inline void TransformPoint(
    OutputPointType opp, const CoefficientPointerVectorType mu,
    const OffsetValueType * gridOffsetTable,
    const double * weights1D )
  {
    for( unsigned int j = 0; j < OutputDimension; ++j )
    {
      opp[ j ] = *( mu[ j ] );
    }
  } // end TransformPoint()
 
 
  static inline void GetJacobian(
    ScalarType * & jacobians, const double * weights1D, double value )
  {
    unsigned long offset = 0;
    for( unsigned int j = 0; j < OutputDimension; ++j )
    {
      offset                  = j * BSplineNumberOfIndices * ( OutputDimension + 1 );
      jacobians[ offset ] = value;
    }
    ++jacobians;
  } // end GetJacobian()
 
 
  static inline void EvaluateJacobianWithImageGradientProduct(
    ScalarType * & imageJacobian, const InternalFloatType * movingImageGradient,
    const double * weights1D, double value )
  {
    for( unsigned int j = 0; j < OutputDimension; ++j )
    {
      *( imageJacobian + j * BSplineNumberOfIndices ) = value * movingImageGradient[ j ];
    }
    ++imageJacobian;
  } // end EvaluateJacobianWithImageGradientProduct()
 
 
  static inline void ComputeNonZeroJacobianIndices(
    unsigned long * & nzji,
    const unsigned long parametersPerDim,
    unsigned long currentIndex,
    const OffsetValueType * gridOffsetTable )
  {
    for( unsigned int j = 0; j < OutputDimension; ++j )
    {
      nzji[ j * BSplineNumberOfIndices ] = currentIndex + j * parametersPerDim;
    }
    ++nzji;
  } // end ComputeNonZeroJacobianIndices()
 
 
  static inline void GetSpatialJacobian(
    InternalFloatType * sj,
    const CoefficientPointerVectorType mu,
    const OffsetValueType * gridOffsetTable,
    const double * weights1D,                    // normal B-spline weights
    const double * derivativeWeights1D )         // 1st derivative of B-spline
  {
    for( unsigned int j = 0; j < OutputDimension; ++j )
    {
      sj[ j ] = *( mu[ j ] );
    }
  } // end GetSpatialJacobian()
 
 
  static inline void GetSpatialHessian(
    InternalFloatType * sh,
    const CoefficientPointerVectorType mu,
    const OffsetValueType * gridOffsetTable,
    const double * weights1D,                   // normal B-spline weights
    const double * derivativeWeights1D,         // 1st derivative of B-spline
    const double * hessianWeights1D )           // 2nd derivative of B-spline
  {
    for( unsigned int j = 0; j < OutputDimension; ++j )
    {
      sh[ j ] = *( mu[ j ] );
    }
  } // end GetSpatialHessian()
 
 
  static inline void GetJacobianOfSpatialJacobian(
    InternalFloatType * & jsj_out,
    const double * weights1D,                   // normal B-spline weights
    const double * derivativeWeights1D,         // 1st derivative of B-spline
    const double * directionCosines,
    InternalFloatType * jsj )
  {
    for( unsigned int j = 0; j < OutputDimension; ++j )
    {
      jsj_out[ j ] = jsj[ OutputDimension ] * directionCosines[ j ];
      for( unsigned int k = 1; k < OutputDimension; ++k )
      {
        jsj_out[ k ] += jsj[ OutputDimension - k ] * directionCosines[ k * OutputDimension + j ];
      }
    }
 
    unsigned int offset = 0;
    for( unsigned int i = 0; i < OutputDimension; ++i )
    {
      offset = i * ( OutputDimension * ( BSplineNumberOfIndices * OutputDimension + 1 ) );
      for( unsigned int j = 0; j < OutputDimension; ++j )
      {
        jsj_out[ j + offset ] = jsj_out[ j ];
      }
    }
 
    jsj_out += OutputDimension * OutputDimension;
 
  } // end GetJacobianOfSpatialJacobian()
 
 
  static inline void GetJacobianOfSpatialHessian(
    InternalFloatType * & jsh_out,
    const double * weights1D,                   // normal B-spline weights
    const double * derivativeWeights1D,         // 1st derivative of B-spline
    const double * hessianWeights1D,            // 2nd derivative of B-spline
    const double * directionCosines,
    InternalFloatType * jsh )
  {
    double jsh_tmp[ OutputDimension * OutputDimension ];
    double matrixProduct[ OutputDimension * OutputDimension ];
 
    if( OutputDimension == 3 )
    {
      jsh_tmp[ 0 ] = jsh[ 9 ];        jsh_tmp[ 1 ] = jsh[ 8 ];        jsh_tmp[ 2 ] = jsh[ 7 ];
      jsh_tmp[ 3 ] = jsh_tmp[ 1 ];    jsh_tmp[ 4 ] = jsh[ 5 ];        jsh_tmp[ 5 ] = jsh[ 4 ];
      jsh_tmp[ 6 ] = jsh_tmp[ 2 ];    jsh_tmp[ 7 ] = jsh_tmp[ 5 ];    jsh_tmp[ 8 ] = jsh[ 2 ];
    }
    else if( OutputDimension == 2 )
    {
      jsh_tmp[ 0 ] = jsh[ 5 ];        jsh_tmp[ 1 ] = jsh[ 4 ];
      jsh_tmp[ 2 ] = jsh_tmp[ 1 ];    jsh_tmp[ 3 ] = jsh[ 2 ];
    }
    else // the general case
    {
      for( unsigned int j = 0; j < OutputDimension; ++j )
      {
        for( unsigned int i = 0; i <= j; ++i )
        {
          jsh_tmp[ j * OutputDimension + i ] = jsh[ ( OutputDimension - j ) + ( OutputDimension - i ) * ( OutputDimension - i + 1 ) / 2 ];
          if( i != j )
          {
            jsh_tmp[ i * OutputDimension + j ] = jsh_tmp[ j * OutputDimension + i ];
          }
        }
      }
    }
 
    for( unsigned int i = 0; i < OutputDimension; ++i )   // row
    {
      for( unsigned int j = 0; j < OutputDimension; ++j ) // column
      {
        double accum = directionCosines[ i ] * jsh_tmp[ j ];
        for( unsigned int k = 1; k < OutputDimension; ++k )
        {
          accum += directionCosines[ k * OutputDimension + i ] * jsh_tmp[ k * OutputDimension + j ];
        }
        matrixProduct[ i * OutputDimension + j ] = accum;
      }
    }
 
    for( unsigned int i = 0; i < OutputDimension; ++i )   // row
    {
      for( unsigned int j = 0; j < OutputDimension; ++j ) // column
      {
        double accum = matrixProduct[ i * OutputDimension ] * directionCosines[ j ];
        for( unsigned int k = 1; k < OutputDimension; ++k )
        {
          accum += matrixProduct[ i * OutputDimension + k ] * directionCosines[ k * OutputDimension + j ];
        }
        jsh_out[ i * OutputDimension + j ] = accum;
      }
    }
 
    unsigned long offset = 0;
    for( unsigned int i = 0; i < OutputDimension; ++i )
    {
      offset = i * ( OutputDimension * OutputDimension * ( BSplineNumberOfIndices * OutputDimension + 1 ) );
      for( unsigned int j = 0; j < OutputDimension * OutputDimension; ++j )
      {
        jsh_out[ j + offset ] = jsh_out[ j ];
      }
    }
 
    jsh_out += OutputDimension * OutputDimension * OutputDimension;
 
  } // end GetJacobianOfSpatialHessian()
 
 
};
 
 
} // end namespace itk
 
#endif /* __itkRecursiveBSplineTransformImplementation_h */