ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkFEMSolver.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 
00019 #ifndef __itkFEMSolver_h
00020 #define __itkFEMSolver_h
00021 
00022 #include "itkProcessObject.h"
00023 #include "itkFEMObject.h"
00024 
00025 #include "itkFEMLinearSystemWrapper.h"
00026 #include "itkFEMLinearSystemWrapperVNL.h"
00027 
00028 #include "itkImage.h"
00029 
00030 namespace itk
00031 {
00032 namespace fem
00033 {
00068 template <unsigned int VDimension = 3>
00069 class ITK_EXPORT Solver : public ProcessObject
00070 {
00071 public:
00072 
00074   typedef Solver                   Self;
00075   typedef ProcessObject            Superclass;
00076   typedef SmartPointer<Self>       Pointer;
00077   typedef SmartPointer<const Self> ConstPointer;
00078 
00080   itkNewMacro(Self);
00081 
00083   itkTypeMacro(Solver, ProcessObject);
00084 
00085   itkStaticConstMacro(FEMDimension, unsigned int, VDimension);
00086   itkStaticConstMacro(MaxDimensions, unsigned int, 3);
00087 
00089   typedef typename itk::fem::FEMObject<VDimension> FEMObjectType;
00090   typedef typename FEMObjectType::Pointer          FEMObjectPointer;
00091   typedef typename FEMObjectType::ConstPointer     FEMObjectConstPointer;
00092   typedef typename DataObject::Pointer             DataObjectPointer;
00093 
00095   typedef Element::Float           Float;
00096   typedef Element::VectorType      VectorType;
00097   typedef Element::Node::ArrayType NodeArray;
00098   typedef Element::ArrayType       ElementArray;
00099   typedef Load::ArrayType          LoadArray;
00100   typedef Material::ArrayType      MaterialArray;
00101 
00105   typedef typename itk::Image<Element::ConstPointer, VDimension> InterpolationGridType;
00106   typedef typename InterpolationGridType::Pointer                InterpolationGridPointerType;
00107   typedef typename InterpolationGridType::SizeType               InterpolationGridSizeType;
00108   typedef typename InterpolationGridType::RegionType             InterpolationGridRegionType;
00109   typedef typename InterpolationGridType::PointType              InterpolationGridPointType;
00110   typedef typename InterpolationGridType::SpacingType            InterpolationGridSpacingType;
00111   typedef typename InterpolationGridType::IndexType              InterpolationGridIndexType;
00112   typedef typename InterpolationGridType::DirectionType          InterpolationGridDirectionType;
00113 
00114   /*
00115    * Get/Set the Interpolation Grid Origin
00116    */
00117   itkSetMacro(Origin, InterpolationGridPointType);
00118   itkGetMacro(Origin, InterpolationGridPointType);
00119 
00120   /*
00121    * Get/Set the Interpolation Grid Spacing
00122    */
00123   itkSetMacro(Spacing, InterpolationGridSpacingType);
00124   itkGetMacro(Spacing, InterpolationGridSpacingType);
00125 
00126   /*
00127    * Get/Set the Interpolation Grid Region
00128    */
00129   itkSetMacro(Region, InterpolationGridRegionType);
00130   itkGetMacro(Region, InterpolationGridRegionType);
00131 
00132   /*
00133    * Get/Set the Interpolation Grid Direction
00134    */
00135   itkSetMacro(Direction, InterpolationGridDirectionType);
00136   itkGetMacro(Direction, InterpolationGridDirectionType);
00137 
00139   virtual Float GetTimeStep(void) const;
00140 
00146   virtual void SetTimeStep(Float dt);
00147 
00149   Float GetSolution(unsigned int i, unsigned int which = 0);
00150 
00153   using Superclass::SetInput;
00154   virtual void SetInput( FEMObjectType *fem);
00155 
00156   virtual void SetInput( unsigned int, FEMObjectType * fem);
00157 
00158   FEMObjectType * GetInput(void);
00159 
00160   FEMObjectType * GetInput(unsigned int idx);
00161 
00170   const Element * GetElementAtPoint(const VectorType & pt) const;
00171 
00173   Float GetDeformationEnergy(unsigned int SolutionIndex = 0);
00174 
00189   void SetLinearSystemWrapper(LinearSystemWrapper::Pointer ls);
00190 
00196   LinearSystemWrapper::Pointer GetLinearSystemWrapper()
00197   {
00198     return m_ls;
00199   }
00200 
00215   void InitializeInterpolationGrid(const InterpolationGridSizeType & size, const InterpolationGridPointType & bb1,
00216                                    const InterpolationGridPointType & bb2);
00217 
00219   void InitializeInterpolationGrid(const InterpolationGridSizeType & size)
00220   {
00221     InterpolationGridPointType bb1;
00223 
00224     bb1.Fill(0.0);
00225 
00226     InterpolationGridPointType bb2;
00227     for( unsigned int i = 0; i < FEMDimension; i++ )
00228     {
00229       bb2[i] = size[i] - 1.0;
00230     }
00231     InitializeInterpolationGrid(size, bb1, bb2);
00232   }
00233 
00237   void InitializeInterpolationGrid(const InterpolationGridRegionType & region,
00238                                    const InterpolationGridPointType & origin,
00239                                    const InterpolationGridSpacingType & spacing,
00240                                    const InterpolationGridDirectionType & direction);
00241 
00252   const InterpolationGridType * GetInterpolationGrid(void) const
00253   {
00254     return m_InterpolationGrid.GetPointer();
00255   }
00256 
00257 
00260   typedef ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType;
00261   using Superclass::MakeOutput;
00262   virtual DataObjectPointer MakeOutput(DataObjectPointerArraySizeType);
00263 
00277   FEMObjectType * GetOutput(void);
00278 
00279   FEMObjectType * GetOutput(unsigned int idx);
00280 
00281 protected:
00282   Solver();
00283   virtual ~Solver();
00284   void PrintSelf(std::ostream& os, Indent indent) const;
00285 
00288   void  GenerateData();
00289 
00290 
00299   // void GenerateGFN(void);
00300 
00304   void AssembleK();
00305 
00312   virtual void InitializeMatrixForAssembly(unsigned int N);
00313 
00319   virtual void FinalizeMatrixAfterAssembly(void)
00320   {
00321     // Apply the boundary conditions to the K matrix
00322     this->ApplyBC();
00323   }
00324 
00331   virtual void AssembleElementMatrix(Element::Pointer e);
00332 
00340   virtual void AssembleLandmarkContribution(Element::ConstPointer e, float);
00341 
00355   void ApplyBC(int dim = 0, unsigned int matrix = 0);
00356 
00364   void AssembleF(int dim = 0);
00365 
00369   void DecomposeK(void);
00370 
00374   virtual void RunSolver(void);
00375 
00380   void UpdateDisplacements(void);
00381 
00385   void FillInterpolationGrid(void);
00386 
00391   virtual void InitializeLinearSystemWrapper(void);
00392 
00396   unsigned int m_NGFN;
00397 
00402   unsigned int m_NMFC;
00403 
00405   LinearSystemWrapper::Pointer m_ls;
00406 
00410   LinearSystemWrapperVNL m_lsVNL;
00411 
00417   InterpolationGridPointerType m_InterpolationGrid;
00418 
00419   FEMObjectPointer m_FEMObject;
00420 private:
00421   Solver(const Self &);         // purposely not implemented
00422   void operator=(const Self &); // purposely not implemented
00423 
00424   /*
00425    * Properties of the Interpolation Grid
00426    */
00427   InterpolationGridRegionType       m_Region;
00428   InterpolationGridPointType        m_Origin;
00429   InterpolationGridSpacingType      m_Spacing;
00430   InterpolationGridDirectionType    m_Direction;
00431 
00432 };
00433 }  // end namespace fem
00434 }  // end namespace itk
00435 
00436 #ifndef ITK_MANUAL_INSTANTIATION
00437 #include "itkFEMSolver.hxx"
00438 #endif
00439 
00440 #endif // #ifndef __itkFEMSolver_h
00441