NegJacobianConstraint.h

/*
 * Medical Image Registration ToolKit (MIRTK)
 *
 * Copyright 2017 Imperial College London
 * Copyright 2017 Andreas Schuh
 *
 * 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
 *
 * 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 MIRTK_NegJacobianConstraint_H
#define MIRTK_NegJacobianConstraint_H

#include "mirtk/JacobianConstraint.h"


namespace mirtk {


/**
 * Penalize small and negative Jacobian determinant values of parameterization
 *
 * This constraint prevents folding of the transformation parameterization,
 * i.e., either of the control point displacements or velocities. It preserves
 * topology in case of a classic FFD model and is in this case identical to
 * the TopologyPreservationConstraint. For the SVFFD model, it penalizes
 * small and negative Jacobian determinant values of the velocity field,
 * not the displacement field generated by this stationary velocity field.
 *
 * This Jacobian based transformation constraint is based on the original
 * IRTK implementation, but not identical. When the Jacobian determinant
 * value is smaller than a small value, _Epsilon, the penalty function is
 * a linear function instead, which is a continuation of the penalty function
 * used for determinant values in the range _Epsilon < det(J) < _Gamma.
 * The slope of this linear function corresponds to the derivative of the
 * penalty at _Epsilon. The problem with the IRTK implementation is that the
 * penalty decreases again for negative values before increasing slowly
 * quadratically. It is further undefined for zero determinant values.
 * This is because the penalty function adopted from Edwards et al. (1998)
 * is only defined for positive values.
 *
 * Edwards et al. (1998). A three-component deformation model for image-guided surgery.
 * Medical Image Analysis, 2(4), 355–367. https://doi.org/10.1016/S1361-8415(98)80016-9
 *
 * Rueckert et al. (2006). Diffeomorphic registration using B-splines. MICCAI, 702–709.
 */
class NegJacobianConstraint : public JacobianConstraint
{
  mirtkEnergyTermMacro(LogJacobianConstraint, EM_NegDetJac);

  /// Small value below which Jacobian determinant is penalised linearly
  /// with increasing smaller (i.e., negative) value
  ///
  /// The penalty implemented in IRTK can be used with _Epsilon set to a very
  /// negative value, i.e., negative with great magnitude, -inf, or NaN.
  mirtkPublicAttributeMacro(double, Epsilon);

  /// Jacobian determinant threshold below which penalty is non-zero
  ///
  /// The _Gamma value that was hard coded in IRTK is 0.3. A higher value
  /// is used by this class, however, because it allows the panlty to earlier
  /// prevent "shocks" in the deformation before these impose a too strong
  /// and sudden penalty.
  mirtkPublicAttributeMacro(double, Gamma);

  /// Power of Jacobian determinant used for penalty term
  mirtkPublicAttributeMacro(int, Power);

public:

  /// Constructor
  NegJacobianConstraint(const char * = "", bool = true);

  /// Destructor
  virtual ~NegJacobianConstraint();

protected:

  /// Set parameter value from string
  virtual bool SetWithPrefix(const char *, const char *);

  /// Set parameter value from string
  virtual bool SetWithoutPrefix(const char *, const char *);

public:

  // Import other overloads
  using JacobianConstraint::Parameter;

  /// Get parameter key/value as string map
  virtual ParameterList Parameter() const;

  // ---------------------------------------------------------------------------
  // Penalty

public:

  /// Evaluate penalty at control point location given Jacobian determinant value
  virtual double Penalty(double det) const
  {
    if (det >= _Gamma) return 0.;
    double a = pow(_Gamma, _Power);
    if (det < _Epsilon) {
      double b = pow(_Epsilon, _Power);
      double m = _Power * ((pow(_Epsilon, _Power - 1) / a) - (a / pow(_Epsilon, _Power + 1)));
      double t = a / b + b / a - 2. - m * _Epsilon;
      return m * det + t;
    } else {
      double b = pow(det, _Power);
      return a / b + b / a - 2.;
    }
  }

  /// Evaluate derivative of penalty w.r.t. Jacobian determinant value
  virtual double DerivativeWrtJacobianDet(double det) const
  {
    if (det >= _Gamma) return 0.;
    if (det < _Epsilon) det = _Epsilon;
    double a = pow(_Gamma, _Power);
    return _Power * ((pow(det, _Power - 1) / a) - (a / pow(det, _Power + 1)));
  }

};


} // namespace mirtk

#endif // MIRTK_NegJacobianConstraint_H