ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkRegionBasedLevelSetFunction.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 #ifndef __itkRegionBasedLevelSetFunction_h
00019 #define __itkRegionBasedLevelSetFunction_h
00020 
00021 #include "itkFiniteDifferenceFunction.h"
00022 #include "itkRegularizedHeavisideStepFunction.h"
00023 #include "vnl/vnl_matrix_fixed.h"
00024 
00025 namespace itk
00026 {
00064 template< class TInput,   // LevelSetImageType
00065           class TFeature, // FeatureImageType
00066           class TSharedData >
00067 class ITK_EXPORT RegionBasedLevelSetFunction:public
00068   FiniteDifferenceFunction< TInput >
00069 {
00070 public:
00072   typedef RegionBasedLevelSetFunction        Self;
00073   typedef FiniteDifferenceFunction< TInput > Superclass;
00074   typedef SmartPointer< Self >               Pointer;
00075   typedef SmartPointer< const Self >         ConstPointer;
00076 
00077   itkStaticConstMacro(ImageDimension, unsigned int, Superclass::ImageDimension);
00078 
00079   // itkNewMacro() is purposely not provided since this is an abstract class.
00080 
00082   itkTypeMacro(RegionBasedLevelSetFunction, FiniteDifferenceFunction);
00083 
00085   typedef double                                      TimeStepType;
00086   typedef typename Superclass::ImageType              ImageType;
00087   typedef typename Superclass::PixelType              PixelType;
00088   typedef PixelType                                   ScalarValueType;
00089   typedef typename Superclass::RadiusType             RadiusType;
00090   typedef typename Superclass::NeighborhoodType       NeighborhoodType;
00091   typedef typename Superclass::NeighborhoodScalesType NeighborhoodScalesType;
00092   typedef typename Superclass::FloatOffsetType        FloatOffsetType;
00093   typedef FixedArray< ScalarValueType, itkGetStaticConstMacro(ImageDimension) >
00094   VectorType;
00095 
00096   /* This structure is derived from LevelSetFunction and stores intermediate
00097   values for computing time step sizes */
00098   struct GlobalDataStruct {
00099     GlobalDataStruct()
00100     {
00101       ScalarValueType null_value = NumericTraits< ScalarValueType >::Zero;
00102 
00103       m_MaxCurvatureChange   = null_value;
00104       m_MaxAdvectionChange   = null_value;
00105       m_MaxGlobalChange      = null_value;
00106     }
00107 
00108     ~GlobalDataStruct() {}
00109 
00110     vnl_matrix_fixed< ScalarValueType,
00111                       itkGetStaticConstMacro(ImageDimension),
00112                       itkGetStaticConstMacro(ImageDimension) > m_dxy;
00113 
00114     ScalarValueType m_dx[itkGetStaticConstMacro(ImageDimension)];
00115 
00116     ScalarValueType m_dx_forward[itkGetStaticConstMacro(ImageDimension)];
00117     ScalarValueType m_dx_backward[itkGetStaticConstMacro(ImageDimension)];
00118 
00119     ScalarValueType m_GradMagSqr;
00120     ScalarValueType m_GradMag;
00121 
00122     ScalarValueType m_MaxCurvatureChange;
00123     ScalarValueType m_MaxAdvectionChange;
00124     ScalarValueType m_MaxGlobalChange;
00125   };
00126 
00127   typedef TInput                                  InputImageType;
00128   typedef typename InputImageType::ConstPointer   InputImageConstPointer;
00129   typedef typename InputImageType::Pointer        InputImagePointer;
00130   typedef typename InputImageType::PixelType      InputPixelType;
00131   typedef typename InputImageType::IndexType      InputIndexType;
00132   typedef typename InputImageType::IndexValueType InputIndexValueType;
00133   typedef typename InputImageType::SizeType       InputSizeType;
00134   typedef typename InputImageType::SizeValueType  InputSizeValueType;
00135   typedef typename InputImageType::RegionType     InputRegionType;
00136   typedef typename InputImageType::PointType      InputPointType;
00137 
00138   typedef TFeature                                FeatureImageType;
00139   typedef typename FeatureImageType::ConstPointer FeatureImageConstPointer;
00140   typedef typename FeatureImageType::PixelType    FeaturePixelType;
00141   typedef typename FeatureImageType::IndexType    FeatureIndexType;
00142   typedef typename FeatureImageType::SpacingType  FeatureSpacingType;
00143   typedef typename FeatureImageType::OffsetType   FeatureOffsetType;
00144 
00145   typedef TSharedData                      SharedDataType;
00146   typedef typename SharedDataType::Pointer SharedDataPointer;
00147 
00148   typedef HeavisideStepFunctionBase< InputPixelType, InputPixelType > HeavisideFunctionType;
00149   typedef typename HeavisideFunctionType::ConstPointer                HeavisideFunctionConstPointer;
00150 
00151   void SetDomainFunction(const HeavisideFunctionType *f)
00152   {
00153     this->m_DomainFunction = f;
00154   }
00155 
00156   virtual void Initialize(const RadiusType & r)
00157   {
00158     this->SetRadius(r);
00159 
00160     // Dummy neighborhood.
00161     NeighborhoodType it;
00162     it.SetRadius(r);
00163 
00164     // Find the center index of the neighborhood.
00165     m_Center =  it.Size() / 2;
00166 
00167     // Get the stride length for each axis.
00168     for ( unsigned int i = 0; i < ImageDimension; i++ )
00169       {
00170       m_xStride[i] = it.GetStride(i);
00171       }
00172   }
00173 
00174 #if !defined( CABLE_CONFIGURATION )
00175   void SetSharedData(SharedDataPointer sharedDataIn)
00176   {
00177     this->m_SharedData = sharedDataIn;
00178   }
00179 #endif
00180 
00181   void UpdateSharedData(bool forceUpdate);
00182 
00183   void * GetGlobalDataPointer() const
00184   {
00185     return new GlobalDataStruct;
00186   }
00187 
00188   TimeStepType ComputeGlobalTimeStep(void *GlobalData) const;
00189 
00191   virtual PixelType ComputeUpdate( const NeighborhoodType & neighborhood,
00192                                    void *globalData, const FloatOffsetType & = FloatOffsetType(0.0) );
00193 
00194   void SetInitialImage(InputImageType *f)
00195   {
00196     m_InitialImage = f;
00197   }
00198 
00199   virtual const FeatureImageType * GetFeatureImage() const
00200   { return m_FeatureImage.GetPointer(); }
00201   virtual void SetFeatureImage(const FeatureImageType *f)
00202   {
00203     m_FeatureImage = f;
00204 
00205     FeatureSpacingType spacing = m_FeatureImage->GetSpacing();
00206     for ( unsigned int i = 0; i < ImageDimension; i++ )
00207       {
00208       this->m_InvSpacing[i] = 1 / spacing[i];
00209       }
00210   }
00211 
00213   virtual VectorType AdvectionField(const NeighborhoodType &,
00214                                     const FloatOffsetType &, GlobalDataStruct * = 0)  const
00215   { return this->m_ZeroVectorConstant; }
00216 
00218   void SetAreaWeight(const ScalarValueType & nu)
00219   { this->m_AreaWeight = nu; }
00220   ScalarValueType GetAreaWeight() const
00221   { return this->m_AreaWeight; }
00223 
00225   void SetLambda1(const ScalarValueType & lambda1)
00226   { this->m_Lambda1 = lambda1; }
00227   ScalarValueType GetLambda1() const
00228   { return this->m_Lambda1; }
00230 
00232   void SetLambda2(const ScalarValueType & lambda2)
00233   { this->m_Lambda2 = lambda2; }
00234   ScalarValueType GetLambda2() const
00235   { return this->m_Lambda2; }
00237 
00239   void SetOverlapPenaltyWeight(const ScalarValueType & gamma)
00240   { this->m_OverlapPenaltyWeight = gamma; }
00241   ScalarValueType GetOverlapPenaltyWeight() const
00242   { return this->m_OverlapPenaltyWeight; }
00244 
00246   virtual void SetCurvatureWeight(const ScalarValueType c)
00247   { m_CurvatureWeight = c; }
00248   ScalarValueType GetCurvatureWeight() const
00249   { return m_CurvatureWeight; }
00251 
00252   void SetAdvectionWeight(const ScalarValueType & iA)
00253   { this->m_AdvectionWeight = iA; }
00254   ScalarValueType GetAdvectionWeight() const
00255   { return this->m_AdvectionWeight; }
00256 
00258   void SetReinitializationSmoothingWeight(const ScalarValueType c)
00259   { m_ReinitializationSmoothingWeight = c; }
00260   ScalarValueType GetReinitializationSmoothingWeight() const
00261   { return m_ReinitializationSmoothingWeight; }
00263 
00265   void SetVolumeMatchingWeight(const ScalarValueType & tau)
00266   { this->m_VolumeMatchingWeight = tau; }
00267   ScalarValueType GetVolumeMatchingWeight() const
00268   { return this->m_VolumeMatchingWeight; }
00270 
00272   void SetVolume(const ScalarValueType & volume)
00273   { this->m_Volume = volume; }
00274   ScalarValueType GetVolume() const
00275   { return this->m_Volume; }
00277 
00279   void SetFunctionId(const unsigned int & iFid)
00280   { this->m_FunctionId = iFid; }
00281 
00282   virtual void ReleaseGlobalDataPointer(void *GlobalData) const
00283   { delete (GlobalDataStruct *)GlobalData; }
00284 
00285   virtual ScalarValueType ComputeCurvature(const NeighborhoodType &,
00286                                            const FloatOffsetType &, GlobalDataStruct *gd);
00287 
00290   virtual ScalarValueType LaplacianSmoothingSpeed(
00291     const NeighborhoodType &,
00292     const FloatOffsetType &, GlobalDataStruct * = 0) const
00293   { return NumericTraits< ScalarValueType >::One; }
00294 
00297   virtual ScalarValueType CurvatureSpeed(const NeighborhoodType &,
00298                                          const FloatOffsetType &, GlobalDataStruct * = 0
00299                                          ) const
00300   { return NumericTraits< ScalarValueType >::One; }
00301 
00306   virtual void CalculateAdvectionImage() {}
00307 protected:
00308 
00309   RegionBasedLevelSetFunction();
00310   virtual ~RegionBasedLevelSetFunction() {}
00311 
00313   InputImageConstPointer m_InitialImage;
00314 
00316   FeatureImageConstPointer m_FeatureImage;
00317 
00318   SharedDataPointer m_SharedData;
00319 
00320   HeavisideFunctionConstPointer m_DomainFunction;
00321 
00323   ScalarValueType m_AreaWeight;
00324 
00326   ScalarValueType m_Lambda1;
00327 
00329   ScalarValueType m_Lambda2;
00330 
00332   ScalarValueType m_OverlapPenaltyWeight;
00333 
00335   ScalarValueType m_VolumeMatchingWeight;
00336 
00338   ScalarValueType m_Volume;
00339 
00341   ScalarValueType m_CurvatureWeight;
00342 
00343   ScalarValueType m_AdvectionWeight;
00344 
00346   ScalarValueType m_ReinitializationSmoothingWeight;
00347 
00348   unsigned int m_FunctionId;
00349 
00350   std::slice x_slice[itkGetStaticConstMacro(ImageDimension)];
00351   OffsetValueType m_Center;
00352   OffsetValueType m_xStride[itkGetStaticConstMacro(ImageDimension)];
00353   double m_InvSpacing[itkGetStaticConstMacro(ImageDimension)];
00354 
00355   static double m_WaveDT;
00356   static double m_DT;
00357 
00358   void ComputeHImage();
00359 
00362   ScalarValueType ComputeGlobalTerm(
00363     const ScalarValueType & imagePixel,
00364     const InputIndexType & inputIndex);
00365 
00370   virtual ScalarValueType ComputeInternalTerm(const FeaturePixelType & iValue,
00371                                               const FeatureIndexType & iIdx) = 0;
00372 
00376   virtual ScalarValueType ComputeExternalTerm(const FeaturePixelType & iValue,
00377                                               const FeatureIndexType & iIdx) = 0;
00378 
00383   virtual ScalarValueType ComputeOverlapParameters(const FeatureIndexType & featIndex,
00384                                                    ScalarValueType & pr) = 0;
00385 
00390   ScalarValueType ComputeVolumeRegularizationTerm();
00391 
00405   ScalarValueType ComputeLaplacian(GlobalDataStruct *gd);
00406 
00408   void ComputeHessian(const NeighborhoodType & it,
00409                       GlobalDataStruct *globalData);
00410 
00412   virtual void ComputeParameters() = 0;
00413 
00416   virtual void UpdateSharedDataParameters() = 0;
00417 
00418   bool m_UpdateC;
00419 
00422   static VectorType InitializeZeroVectorConstant();
00423 
00425   static VectorType m_ZeroVectorConstant;
00426 private:
00427   RegionBasedLevelSetFunction(const Self &); //purposely not implemented
00428   void operator=(const Self &);              //purposely not implemented
00429 };
00430 } // end namespace itk
00432 
00433 #ifndef ITK_MANUAL_INSTANTIATION
00434 #include "itkRegionBasedLevelSetFunction.hxx"
00435 #endif
00436 
00437 #endif
00438