itkFEMRobustSolver.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  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 itkFEMRobustSolver_h
#define itkFEMRobustSolver_h

#include "itkFEMSolver.h"
#include "itkFEMElementBase.h"
#include "itkFEMMaterialBase.h"
#include "itkFEMLoadBase.h"
#include "itkFEMLoadNoisyLandmark.h"
#include "itkFEMLinearSystemWrapperVNL.h"
#include "itkFEMLinearSystemWrapperItpack.h"

#include <cmath>

namespace itk
{
namespace fem
{

template <unsigned int VDimension = 3>
class RobustSolver : public Solver<VDimension>
{
public:
  typedef RobustSolver              Self;
  typedef Solver<VDimension>        Superclass;
  typedef SmartPointer<Self>        Pointer;
  typedef SmartPointer<const Self>  ConstPointer;

  itkNewMacro(Self);

  itkTypeMacro(RobustSolver, Solver);

  typedef typename Superclass::VectorType                           VectorType;
  typedef typename Superclass::Float                                Float;
  typedef typename Superclass::InterpolationGridType                InterpolationGridType;
  typedef typename Superclass::InterpolationGridPointerType         InterpolationGridPointerType;
  typedef typename Superclass::InterpolationGridSizeType            InterpolationGridSizeType;
  typedef typename Superclass::InterpolationGridRegionType          InterpolationGridRegionType;
  typedef typename Superclass::InterpolationGridPointType           InterpolationGridPointType;
  typedef typename Superclass::InterpolationGridSpacingType         InterpolationGridSpacingType;
  typedef typename Superclass::InterpolationGridIndexType           InterpolationGridIndexType;
  typedef typename InterpolationGridType::DirectionType             InterpolationGridDirectionType;

  itkStaticConstMacro(FEMDimension, unsigned int, VDimension);

  typedef typename Superclass::FEMObjectType    FEMObjectType;

  typedef typename Element::MatrixType                     MatrixType;
  typedef typename FEMObjectType::LoadContainerType        LoadContainerType;
  typedef typename FEMObjectType::NodeContainerType        NodeContainerType;
  typedef typename FEMObjectType::LoadContainerIterator    LoadContainerIterator;

  itkSetMacro(ApproximationSteps, unsigned int);
  itkGetMacro(ApproximationSteps, unsigned int);

  itkSetMacro(OutlierRejectionSteps, unsigned int);
  itkGetMacro(OutlierRejectionSteps, unsigned int);

  itkSetMacro(ToleranceToLargestDisplacement, double);
  itkGetMacro(ToleranceToLargestDisplacement, double);

  itkSetMacro(TradeOffImageMeshEnergy, double);
  itkGetMacro(TradeOffImageMeshEnergy, double);

  itkSetMacro(FractionErrorRejected, double);
  itkGetMacro(FractionErrorRejected, double);

  itkSetMacro(UseInterpolationGrid, bool);
  itkGetMacro(UseInterpolationGrid, bool);

protected:

  RobustSolver();
  ~RobustSolver();

  virtual void GenerateData() ITK_OVERRIDE;

  virtual void RunSolver(void) ITK_OVERRIDE;

  void Initialization();

  void InitializeInterpolationGrid();

  void InitializeLandmarks();

  void AssembleMeshStiffnessMatrix();

  virtual void AssembleElementMatrixWithID(const Element::Pointer & e, unsigned int matrixIndex);

  void AssembleLandmarkStiffnessMatrix();

  void AssembleGlobalMatrixFromLandmarksAndMeshMatrices();

  void AssembleF();

  void IncrementalSolverWithOutlierRejection();

  void IncrementalSolverWithoutOutlierRejection();

  void SolveSystem();

  void ComputeLandmarkSimulatedDisplacementAndWeightedError();

  void ComputeLandmarkTensor();

  float GetLandmarkTensorPonderation() const;

  void NthElementWRTDisplacementError(unsigned int numberOfRejectedBlocksPerStep);

  void UnselectLandmarks(unsigned int numberOfRejectedBlocksPerStep);

  void RemoveUnselectedLandmarkContributionInPointStiffnessMatrix();

  void DeleteFromLandmarkBeginning(unsigned int numberOfRejectedLandmarksPerStep);

  void DeleteLandmarksOutOfMesh();

  void RescaleLandmarkStiffnessMatrix(double oldPointTensorPonderation);

  void CalculateExternalForces();

  void AddExternalForcesToSetMeshZeroEnergy();

private:

  RobustSolver(const Self &) ITK_DELETE_FUNCTION;
  void operator=(const Self &) ITK_DELETE_FUNCTION;

  unsigned int m_OutlierRejectionSteps;

  unsigned int m_ApproximationSteps;


  typedef unsigned int FEMIndexType;

  FEMIndexType m_ForceIndex;
  FEMIndexType m_LandmarkForceIndex;
  FEMIndexType m_ExternalForceIndex;
  FEMIndexType m_SolutionIndex;
  FEMIndexType m_MeshStiffnessMatrixIndex;
  FEMIndexType m_LandmarkStiffnessMatrixIndex;
  FEMIndexType m_StiffnessMatrixIndex;

  double    m_TradeOffImageMeshEnergy;

  double    m_ToleranceToLargestDisplacement;
  double    m_ConjugateGradientPrecision;
  double    m_FractionErrorRejected;

  bool m_UseInterpolationGrid;

  LinearSystemWrapperItpack m_Itpack;
};

class CompareLandmarkDisplacementError :
    public std::binary_function<const Load::Pointer &, const Load::Pointer &, bool>
{
public:
    bool operator()(const Load::Pointer & L1 , const Load::Pointer & L2)
    {
    LoadNoisyLandmark * l1 = dynamic_cast<LoadNoisyLandmark*>(L1.GetPointer());
    LoadNoisyLandmark * l2 = dynamic_cast<LoadNoisyLandmark*>(L2.GetPointer());

    return l1->GetErrorNorm() > l2->GetErrorNorm();
    }
};

}  // end namespace fem
}  // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkFEMRobustSolver.hxx"
#endif

#endif