/* * Medical Image Registration ToolKit (MIRTK) * * Copyright 2013-2017 Imperial College London * Copyright 2013-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_GenericRegistrationFilter_H #define MIRTK_GenericRegistrationFilter_H #include "mirtk/RegistrationFilter.h" #include "mirtk/Math.h" #include "mirtk/Point.h" #include "mirtk/Vector3D.h" #include "mirtk/EventDelegate.h" #include "mirtk/BaseImage.h" #include "mirtk/GenericImage.h" #include "mirtk/RegisteredImage.h" #include "mirtk/InterpolationMode.h" #include "mirtk/ExtrapolationMode.h" #include "mirtk/SimilarityMeasure.h" #include "mirtk/PointSetDistanceMeasure.h" #include "mirtk/ConstraintMeasure.h" #include "mirtk/LocalOptimizer.h" #include "mirtk/RegistrationEnergy.h" #include "mirtk/Transformation.h" #include "mirtk/TransformationModel.h" #include "mirtk/RegistrationConfig.h" #if MIRTK_Registration_WITH_PointSet # include "vtkSmartPointer.h" # include "vtkPointSet.h" #else class vtkPointSet; #endif namespace mirtk { class RegisteredPointSet; class RegistrationEnergyParser; class GenericRegistrationLogger; class GenericRegistrationDebugger; class HomogeneousTransformation; class FreeFormTransformation; const int MAX_NO_RESOLUTIONS = 10; /** * Generic registration filter */ class GenericRegistrationFilter : public RegistrationFilter { mirtkObjectMacro(GenericRegistrationFilter); friend class RegistrationEnergyParser; friend class GenericRegistrationLogger; friend class GenericRegistrationDebugger; // --------------------------------------------------------------------------- // Types public: /// Type of resolution pyramid images typedef RegisteredImage::InputImageType ResampledImageType; /// List type storing images for one resolution pyramid level typedef Array ResampledImageList; /// Scalar type of resolution pyramid images typedef ResampledImageType::VoxelType VoxelType; /// Type of cached displacement field typedef RegisteredImage::DisplacementImageType DisplacementImageType; /// Structure storing information about transformation instance struct TransformationInfo { double _Exponent; TransformationInfo(double e = .0) : _Exponent(e) {} bool operator ==(const TransformationInfo &other) const { return fequal(_Exponent, other._Exponent, 1e-3); } bool operator !=(const TransformationInfo &other) const { return !(*this == other); } static TransformationInfo Identity() { return TransformationInfo(.0); } static TransformationInfo Full() { return TransformationInfo(1.0); } static TransformationInfo Inverse() { return TransformationInfo(-1.0); } bool IsForwardTransformation() const { return _Exponent > .0; } bool IsBackwardTransformation() const { return _Exponent < .0; } bool IsIdentity() const { return fequal(_Exponent, .0, 1e-3); } operator bool() const { return !IsIdentity(); } }; /// Structure storing information of a used image similarity term parsed /// from the registration energy function string struct ImageSimilarityInfo { SimilarityMeasure _Measure; ///< Type of similarity measure bool _DefaultSign; ///< Whether to use default sign of similarity string _Name; ///< Name of similarity term double _Weight; ///< Weight of similarity term int _TargetIndex; ///< Index of target image TransformationInfo _TargetTransformation; ///< Target transformation identifier int _SourceIndex; ///< Index of source image TransformationInfo _SourceTransformation; ///< Source transformation identifier bool IsSymmetric() const { return _TargetTransformation && _SourceTransformation; } }; /// Structure storing information of a used point set distance term parsed /// from the registration energy function string struct PointSetDistanceInfo { PointSetDistanceMeasure _Measure; ///< Measure of polydata distance bool _DefaultSign; ///< Whether to use default sign of distance measure string _Name; ///< Name of polydata distance term double _Weight; ///< Weight of polydata distance term int _TargetIndex; ///< Index of target data set TransformationInfo _TargetTransformation; ///< Target transformation identifier int _SourceIndex; ///< Index of source data set TransformationInfo _SourceTransformation; ///< Source transformation identifier bool IsSymmetric() const { return _TargetTransformation && _SourceTransformation; } }; /// Structure storing information of a used point set constraint term parsed /// from the registration energy function string struct PointSetConstraintInfo { EnergyMeasure _Measure; ///< Type of internal forces string _Name; ///< Name of constraint double _Weight; ///< Weight of constraint int _PointSetIndex; ///< Index of input point set object int _RefImageIndex; int _RefPointSetIndex; TransformationInfo _Transformation; ///< Point set transformation identifier }; /// Structure storing information of a used constraint term parsed from the /// registration energy function string struct ConstraintInfo { ConstraintMeasure _Measure; ///< Type of constraint string _Name; ///< Name of constraint double _Weight; ///< Weight of constraint }; /// Structure storing information about cached displacements struct DisplacementInfo { int _DispIndex; ///< Index of cached displacement field double _InputTime; ///< Time of untransformed data ImageAttributes _Domain; ///< Domain on which it is defined const Transformation *_Transformation; ///< Corresponding transformation instance }; /// Structure storing information about transformed output point set struct PointSetOutputInfo { int _InputIndex; bool _InitialUpdate; TransformationInfo _Transformation; }; // --------------------------------------------------------------------------- // Attributes /// Number of resolution levels mirtkPublicAttributeMacro(int, NumberOfLevels); /// Level at which to stop multi-resolution optimization mirtkPublicAttributeMacro(int, FinalLevel); /// Multi-level transformation mode mirtkPublicAttributeMacro(MFFDMode, MultiLevelMode); /// Transformation model mirtkPublicAttributeMacro(Array, TransformationModel); /// Default image interpolation mode mirtkPublicAttributeMacro(enum InterpolationMode, DefaultInterpolationMode); /// Default image extrapolation mode mirtkPublicAttributeMacro(enum ExtrapolationMode, DefaultExtrapolationMode); /// Rescale input image intensities from [imin, imax] to [0, omax] /// /// Here [imin, imax] is the intensity range of the respective input image /// excluding the background (if specified) and omax is the maximum intensity /// of the rescaled intensities. This maximum value influences the relative /// weighting of image similarity gradient vs. the gradient of constraint terms /// such as the bending energy. /// /// This normalization ensures that the magnitude of image derivatives used /// for energy gradient computation have comparable magnitude in order to /// prevent one image to have stronger influence than another. This is more /// important for symmetric, inverse consistent, or multi-modal settings. /// /// By default, this parameter is set to inf and no rescaling is done. mirtkPublicAttributeMacro(double, MaxRescaledIntensity); /// Whether to precompute image derivatives or compute them on the fly mirtkPublicAttributeMacro(bool, PrecomputeDerivatives); /// Default similarity measure mirtkPublicAttributeMacro(enum SimilarityMeasure, SimilarityMeasure); /// Default polydata distance measure mirtkPublicAttributeMacro(enum PointSetDistanceMeasure, PointSetDistanceMeasure); /// Optimization method mirtkPublicAttributeMacro(enum OptimizationMethod, OptimizationMethod); /// Normalize weights of energy function terms mirtkPublicAttributeMacro(bool, NormalizeWeights); /// Initial guess of optimal transformation mirtkPublicAggregateMacro(const Transformation, InitialGuess); /// Target transformation to be approximated by output transformation mirtkPublicAggregateMacro(const Transformation, TargetTransformation); /// Name of target transformation MSDE term mirtkPublicAttributeMacro(string, TargetTransformationErrorName); /// Weight of target transformation MSDE term mirtkPublicAttributeMacro(double, TargetTransformationErrorWeight); /// Whether to merge initial global transformation into (local) transformation mirtkPublicAttributeMacro(bool, MergeGlobalAndLocalTransformation); /// Whether to allow x coordinate transformation mirtkPublicAttributeMacro(bool, RegisterX); /// Whether to allow y coordinate transformation mirtkPublicAttributeMacro(bool, RegisterY); /// Whether to allow z coordinate transformation mirtkPublicAttributeMacro(bool, RegisterZ); /// Enforce Dirichlet boundary condition on FFD transformations /// /// When this option is enabled, the status of control points at the /// boundary of the finite FFD lattice is set to Passive and the parameters /// of these control points set to zero. This is always the case for FFDs /// whose parameters are control point displacements. In case of FFDs /// parameterized by (stationary) velocities, the default extrapolation /// mode is nearest neighbor, however, and a layer of passive CPs with /// constant value is needed if the velocity should be forced to zero /// outside the finite domain on which the velocity field is defined. /// Alternatively, FFD extrapolation mode "Const" can be used. mirtkPublicAttributeMacro(bool, DirichletBoundaryCondition); /// Mask which defines where to evaluate the energy function mirtkPublicAggregateMacro(BinaryImage, Domain); /// Whether to adaptively remesh surfaces before each gradient step mirtkPublicAttributeMacro(bool, AdaptiveRemeshing); protected: /// Common attributes of (untransformed) input target data sets /// /// These attributes are in particular used to initialize the control point /// grid of a free-form deformation such that the grid is large enough to /// be valid for every target point for which the transformation will be /// evaluated. It is therefore computed from the attributes of the input /// data sets rather than the downsampled data. If a _Domain mask is set, /// the attributes of this mask are copied instead. struct ImageAttributes _RegistrationDomain; Array _Centroid; ///< Centroids of images Point _TargetOffset; ///< Target origin offset Point _SourceOffset; ///< Source origin offset Array _Input; ///< Input images Array _Image; ///< Resolution pyramid Array _Mask; ///< Domain on which to evaluate similarity Transformation *_Transformation; ///< Current estimate Array _TransformationInfo; ///< Meta-data of partial transformation Array _TransformationInstance; ///< Partial transformations Array _InterpolationMode; ///< Interpolation mode of each input image Array _ExtrapolationMode; ///< Extrapolation mode of each input image RegistrationEnergy _Energy; ///< Registration energy LocalOptimizer *_Optimizer; ///< Used optimizer enum TransformationModel _CurrentModel; ///< Current transformation model int _CurrentLevel; ///< Current resolution level EventDelegate _EventDelegate; ///< Forwards optimization events to observers string _EnergyFormula; ///< Registration energy formula as string Array _ImageSimilarityInfo; ///< Parsed similarity measure(s) Array _ConstraintInfo; ///< Parsed constraint(s) Array > _Resolution[MAX_NO_RESOLUTIONS]; ///< Image resolution in mm Array _MinEdgeLength[MAX_NO_RESOLUTIONS]; ///< Minimum edge length in mm Array _MaxEdgeLength[MAX_NO_RESOLUTIONS]; ///< Maximum edge length in mm int _UseGaussianResolutionPyramid; ///< Whether resolution levels correspond to a Gaussian pyramid Array _Blurring[MAX_NO_RESOLUTIONS]; ///< Image blurring value double _DefaultBackground; ///< Default background value Array _Background; ///< Image background value bool _DownsampleWithPadding; ///< Whether to take background into account ///< during initialization of the image pyramid bool _CropPadImages; ///< Whether to crop/pad input images int _CropPadFFD; ///< Whether to crop/pad FFD lattice #if MIRTK_Registration_WITH_PointSet Array _PointSetTime; ///< Time point of input points, curves, and/or surfaces Array > _PointSetInput; ///< Input points, curves, and/or surfaces Array > > _PointSet; ///< Remeshed/-sampled points, curves, and/or surfaces Array _PointSetOutput; ///< Output points, curves, and/or surfaces Array _PointSetOutputInfo; ///< Meta-data of output points, curves, and/or surfaces Array _PointSetDistanceInfo; ///< Parsed point set distance measure(s) Array _PointSetConstraintInfo; ///< Parsed point set constraint measure(s) #else // MIRTK_Registration_WITH_PointSet Array _PointSetTime; ///< Unused Array _PointSetInput; ///< Unused Array > _PointSet; ///< Unused Array _PointSetOutput; ///< Unused Array _PointSetOutputInfo; ///< Unused Array _PointSetDistanceInfo; ///< Unused Array _PointSetConstraintInfo; ///< Unused #endif // MIRTK_Registration_WITH_PointSet int _Centering [MAX_NO_RESOLUTIONS]; ///< Whether to center foreground (if applicable) double _MinControlPointSpacing[MAX_NO_RESOLUTIONS][4]; ///< Control point spacing for FFDs double _MaxControlPointSpacing[MAX_NO_RESOLUTIONS][4]; ///< Control point spacing for FFDs bool _Subdivide [MAX_NO_RESOLUTIONS][4]; ///< Whether to subdivide FFD /// Parameters not considered by the registration filter itself /// /// These parameters are passed on to the sub-modules of the registration /// filter such as the image similarity measure(s), the registration /// constraint(s), and the optimizer. This way, the registration filter does /// not need to know which parameters its sub-modules accept and is better /// decoupled from the implementation of the respective image similarities, /// constraints, and optimizers. ParameterList _Parameter[MAX_NO_RESOLUTIONS]; /// Delegate of pre-update function RegistrationEnergy::PreUpdateFunctionType _PreUpdateDelegate; /// Info of cached displacement fields Array _DisplacementInfo; Array _DisplacementField; // --------------------------------------------------------------------------- // Access to managed data objects and their attributes protected: /// Determine number of temporal frames and temporal sampling attributes virtual int NumberOfFrames(double * = NULL, double * = NULL, double * = NULL) const; /// Get attributes of n-th input image at specified resolution level virtual struct ImageAttributes ImageAttributes(int, int = -1) const; /// Common attributes of (untransformed) input target data sets at specified resolution level virtual struct ImageAttributes RegistrationDomain(int = -1) const; /// Average resolution of (untransformed) target data sets at given level virtual Vector3D AverageOutputResolution(int = -1) const; /// Get (partial/inverse) output transformation virtual Transformation *OutputTransformation(TransformationInfo); /// Initialize output image corresponding to registered (transformed) image virtual void SetInputOf(RegisteredImage *, const struct ImageAttributes &, int, TransformationInfo); /// Get (new) registered output point set /// /// \note Use only when MIRTK_Registration_WITH_PointSet is 1. RegisteredPointSet *OutputPointSet(int, double, TransformationInfo); // --------------------------------------------------------------------------- // Construction/Destruction private: /// Copy constructor /// \note Intentionally not implemented GenericRegistrationFilter(const GenericRegistrationFilter &); /// Assignment operator /// \note Intentionally not implemented void operator =(const GenericRegistrationFilter &); public: /// Constructor GenericRegistrationFilter(); /// Reset filter settings, but keep input virtual void Reset(); /// Reset filter settings and input virtual void Clear(); /// Destructor virtual ~GenericRegistrationFilter(); // --------------------------------------------------------------------------- // Input images /// Set input images of the registration filter void Input(const BaseImage *, const BaseImage *); /// Set input images of the registration filter void Input(int, const BaseImage **); /// Set input images of the registration filter template void Input(int, const GenericImage **); /// Add filter input image /// /// For registration filters which support multiple target/source images /// such as multiple channels and/or frames of a temporal image sequence. /// How the transformation is applied to which inputs depends on the /// respective transformation model and image registration method. void AddInput(const BaseImage *); /// Number of input images int NumberOfImages() const; /// Number of required input images /// \note Only valid after ParseEnergyFormula has been called! int NumberOfRequiredImages() const; /// Determine whether the specified input image either remains untransformed, /// or is transformed by the inverse transformation or a part of it such /// as in case of an inverse consistent and/or symmetric energy function. bool IsTargetImage(int) const; /// Determine whether the specified input image will be transformed by the /// forward transformation or part of it such as in case of a symmetric energy. bool IsSourceImage(int) const; /// Determine whether the specified input image remains untransformed. bool IsFixedImage(int) const; /// Determine whether the specified input image will be transformed. bool IsMovingImage(int) const; /// Set common/default interpolation mode for all input images void InterpolationMode(enum InterpolationMode); /// Set interpolation mode of n-th input image void InterpolationMode(int, enum InterpolationMode); /// Get interpolation mode of n-th input image or default mode, respectively enum InterpolationMode InterpolationMode(int = -1) const; /// Set common/default extrapolation mode for all input images void ExtrapolationMode(enum ExtrapolationMode); /// Set extrapolation mode of n-th input image void ExtrapolationMode(int, enum ExtrapolationMode); /// Get extrapolation mode of n-th input image or default mode, respectively enum ExtrapolationMode ExtrapolationMode(int = -1) const; /// Get background value of n-th input image (after registration done) double BackgroundValue(int) const; // --------------------------------------------------------------------------- // Input simplicial complexes (points, curves, surfaces, tetrahedral meshes) /// Set input point sets of the registration filter /// /// \note Use only when MIRTK_Registration_WITH_PointSet is 1. void Input(vtkPointSet *, vtkPointSet *, double = .0, double = 1.0); /// Set input point set of the registration filter /// /// \note Use only when MIRTK_Registration_WITH_PointSet is 1. void Input(int, vtkPointSet **, double * = NULL); /// Add filter input point set /// /// \note Use only when MIRTK_Registration_WITH_PointSet is 1. void AddInput(vtkPointSet *, double = .0); /// Number of input point sets int NumberOfPointSets() const; /// Number of required input point sets /// \note Only valid after ParseEnergyFormula has been called! int NumberOfRequiredPointSets() const; /// Determine whether the specified input point set will be transformed by the /// forward transformation or part of it such as in case of a symmetric energy. bool IsTargetPointSet(int) const; /// Determine whether the specified input point set either remains untransformed, /// or is transformed by the inverse transformation or a part of it such /// as in case of an inverse consistent and/or symmetric energy function. bool IsSourcePointSet(int) const; /// Determine whether the specified input point set remains untransformed. bool IsFixedPointSet(int) const; /// Determine whether the specified input point set will be transformed. bool IsMovingPointSet(int) const; // --------------------------------------------------------------------------- // Parameter using RegistrationFilter::Read; using RegistrationFilter::Parameter; /// Set (single) transformation model virtual void TransformationModel(enum TransformationModel); /// Parse registration energy function /// \note Optional to call explicitly as this method is called by GuessParameter. virtual void ParseEnergyFormula(int = -1, int = -1, int = -1); /// Guess proper setting for any yet unset parameter /// \note Optional to call explicitly as this method is called by Run. virtual void GuessParameter(); /// Read registration parameters from input stream virtual bool Read(istream &, bool = false); /// Set named parameter from value as string virtual bool Set(const char *, const char *); /// Set named parameter from value as string virtual bool Set(const char *, const char *, int); /// Get parameters as key/value as string map virtual ParameterList Parameter() const; /// Get parameters as key/value as string map virtual ParameterList Parameter(int) const; /// Write registration parameters to file virtual void Write(const char *) const; // --------------------------------------------------------------------------- // Execution /// Make an initial guess of the (global) output transformation virtual Transformation *MakeInitialGuess(); /// Run the multi-level registration virtual void Run(); // --------------------------------------------------------------------------- // Implementation - the following functions can be overridden in subclasses protected: /// Whether current level is the initial resolution level bool AtInitialLevel() const; /// Whether current level is the final resolution level bool AtFinalLevel() const; /// Run multi-resolution registration virtual void MultiResolutionOptimization(); /// Initialize registration at current resolution virtual void Initialize(); /// Initialize image resolution pyramid virtual void InitializePyramid(); /// Remesh/-sample input point sets virtual void InitializePointSets(); /// Type of output transformation of sub-registration at current resolution virtual enum TransformationType TransformationType(); /// Initialize new transformation instance virtual void InitializeTransformation(); /// Initialize transformation parameters using provided initial guess virtual void ApplyInitialGuess(); /// Initialize status of linear parameters virtual void InitializeStatus(HomogeneousTransformation *); /// Initialize status of FFD parameters virtual void InitializeStatus(FreeFormTransformation *); /// Initialize status of transformation parameters virtual void InitializeStatus(); /// Initialize transformation for sub-registration at current resolution virtual void InitializeOutput(); /// Instantiate new image similarity term for energy function /// \note The individual energy terms are destroyed by the energy function! virtual void AddImageSimilarityTerm(); /// Instantiate new point set distance term for energy function /// \note The individual energy terms are destroyed by the energy function! virtual void AddPointSetDistanceTerm(); /// Instantiate new point set constraint term for energy function /// \note The individual energy terms are destroyed by the energy function! virtual void AddPointSetConstraintTerm(); /// Instantiate new regularization term for energy function /// \note The individual energy terms are destroyed by the energy function! virtual void AddPenaltyTerm(); /// Initialize registration energy of registration at current resolution virtual void InitializeEnergy(); /// Initialize optimizer used to solve registration problem virtual void InitializeOptimizer(); /// Finalize registration at current resolution virtual void Finalize(); /// Callback function called by _Energy->Update(bool) void PreUpdateCallback(bool); }; //////////////////////////////////////////////////////////////////////////////// // Inline/template definitions //////////////////////////////////////////////////////////////////////////////// // ----------------------------------------------------------------------------- template void GenericRegistrationFilter::Input(int num, const GenericImage **image) { _Input.clear(); for (int n = 0; n < num; ++n) AddInput(image[n]); } } // namespace mirtk #endif // MIRTK_GenericRegistrationFilter_H