00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
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
00062 template < class TInput,
00063 class TFeature,
00064 class TSharedData >
00065 class ITK_EXPORT RegionBasedLevelSetFunction: public
00066 FiniteDifferenceFunction< TInput >
00067 {
00068 public:
00070 typedef RegionBasedLevelSetFunction Self;
00071 typedef FiniteDifferenceFunction< TInput > Superclass;
00072 typedef SmartPointer<Self> Pointer;
00073 typedef SmartPointer<const Self> ConstPointer;
00074
00075 itkStaticConstMacro(ImageDimension, unsigned int, Superclass::ImageDimension);
00076
00077
00078
00080 itkTypeMacro( RegionBasedLevelSetFunction, FiniteDifferenceFunction );
00081
00083 typedef double TimeStepType;
00084 typedef typename Superclass::ImageType ImageType;
00085 typedef typename Superclass::PixelType PixelType;
00086 typedef PixelType ScalarValueType;
00087 typedef typename Superclass::RadiusType RadiusType;
00088 typedef typename Superclass::NeighborhoodType NeighborhoodType;
00089 typedef typename Superclass::NeighborhoodScalesType NeighborhoodScalesType;
00090 typedef typename Superclass::FloatOffsetType FloatOffsetType;
00091 typedef FixedArray< ScalarValueType, itkGetStaticConstMacro(ImageDimension) >
00092 VectorType;
00093
00094
00095
00096 struct GlobalDataStruct
00097 {
00098 GlobalDataStruct()
00099 {
00100 ScalarValueType null_value = NumericTraits<ScalarValueType>::Zero;
00101
00102 m_MaxCurvatureChange = null_value;
00103 m_MaxGlobalChange = null_value;
00104 }
00105
00106 ~GlobalDataStruct() {}
00107
00108
00109
00110 ScalarValueType m_MaxCurvatureChange;
00111
00112 vnl_matrix_fixed<ScalarValueType,
00113 itkGetStaticConstMacro(ImageDimension),
00114 itkGetStaticConstMacro(ImageDimension)> m_dxy;
00115
00116 ScalarValueType m_dx[itkGetStaticConstMacro(ImageDimension)];
00117
00118 ScalarValueType m_dx_forward[itkGetStaticConstMacro(ImageDimension)];
00119 ScalarValueType m_dx_backward[itkGetStaticConstMacro(ImageDimension)];
00120
00121 ScalarValueType m_GradMagSqr;
00122 ScalarValueType m_MaxGlobalChange;
00123 };
00124
00125
00126 typedef TInput InputImageType;
00127 typedef typename InputImageType::ConstPointer InputImageConstPointer;
00128 typedef typename InputImageType::Pointer InputImagePointer;
00129 typedef typename InputImageType::PixelType InputPixelType;
00130 typedef typename InputImageType::IndexType InputIndexType;
00131 typedef typename InputImageType::IndexValueType InputIndexValueType;
00132 typedef typename InputImageType::SizeType InputSizeType;
00133 typedef typename InputImageType::SizeValueType InputSizeValueType;
00134 typedef typename InputImageType::RegionType InputRegionType;
00135 typedef typename InputImageType::PointType InputPointType;
00136
00137 typedef TFeature FeatureImageType;
00138 typedef typename FeatureImageType::ConstPointer FeatureImageConstPointer;
00139 typedef typename FeatureImageType::PixelType FeaturePixelType;
00140 typedef typename FeatureImageType::IndexType FeatureIndexType;
00141 typedef typename FeatureImageType::OffsetType FeatureOffsetType;
00142
00143 typedef TSharedData SharedDataType;
00144 typedef typename SharedDataType::Pointer SharedDataPointer;
00145
00146 typedef HeavisideStepFunctionBase< InputPixelType, InputPixelType > HeavisideFunctionType;
00147 typedef typename HeavisideFunctionType::ConstPointer HeavisideFunctionConstPointer;
00148
00149 void SetDomainFunction( const HeavisideFunctionType * f )
00150 {
00151 this->m_DomainFunction = f;
00152 }
00153
00154 virtual void Initialize(const RadiusType &r)
00155 {
00156 this->SetRadius(r);
00157
00158
00159 NeighborhoodType it;
00160 it.SetRadius( r );
00161
00162
00163 m_Center = it.Size() / 2;
00164
00165
00166 for(unsigned int i = 0; i < ImageDimension; i++)
00167 {
00168 m_xStride[i] = it.GetStride(i);
00169 }
00170 }
00171
00172
00173 void SetSharedData( SharedDataPointer sharedDataIn )
00174 {
00175 this->m_SharedData = sharedDataIn;
00176 }
00177
00178 void UpdateSharedData( bool forceUpdate );
00179
00180 void *GetGlobalDataPointer() const
00181 {
00182 return new GlobalDataStruct;
00183 }
00184
00185 TimeStepType ComputeGlobalTimeStep(void *GlobalData) const;
00186
00188 PixelType ComputeUpdate(const NeighborhoodType &neighborhood,
00189 void *globalData, const FloatOffsetType& = FloatOffsetType(0.0));
00190
00191 void SetInitialImage(InputImageType *f)
00192 {
00193 m_InitialImage = f;
00194 }
00195
00196 virtual const FeatureImageType *GetFeatureImage() const
00197 { return m_FeatureImage.GetPointer(); }
00198 virtual void SetFeatureImage(const FeatureImageType *f)
00199 { m_FeatureImage = f; }
00200
00202 void SetAreaWeight( const ScalarValueType& nu)
00203 { this->m_AreaWeight = nu; }
00204 ScalarValueType GetAreaWeight() const
00205 { return this->m_AreaWeight; }
00207
00209 void SetLambda1( const ScalarValueType& lambda1 )
00210 { this->m_Lambda1 = lambda1; }
00211 ScalarValueType GetLambda1() const
00212 { return this->m_Lambda1; }
00214
00216 void SetLambda2( const ScalarValueType& lambda2 )
00217 { this->m_Lambda2 = lambda2; }
00218 ScalarValueType GetLambda2() const
00219 { return this->m_Lambda2; }
00221
00223 void SetOverlapPenaltyWeight( const ScalarValueType& gamma )
00224 { this->m_OverlapPenaltyWeight = gamma; }
00225 ScalarValueType GetOverlapPenaltyWeight() const
00226 { return this->m_OverlapPenaltyWeight; }
00228
00230 virtual void SetCurvatureWeight(const ScalarValueType c)
00231 { m_CurvatureWeight = c; }
00232 ScalarValueType GetCurvatureWeight() const
00233 { return m_CurvatureWeight; }
00235
00237 void SetLaplacianSmoothingWeight(const ScalarValueType c)
00238 { m_LaplacianSmoothingWeight = c; }
00239 ScalarValueType GetLaplacianSmoothingWeight() const
00240 { return m_LaplacianSmoothingWeight; }
00242
00244 void SetVolumeMatchingWeight( const ScalarValueType& tau )
00245 { this->m_VolumeMatchingWeight = tau; }
00246 ScalarValueType GetVolumeMatchingWeight() const
00247 { return this->m_VolumeMatchingWeight; }
00249
00251 void SetVolume( const ScalarValueType& volume )
00252 { this->m_Volume = volume; }
00253 ScalarValueType GetVolume() const
00254 { return this->m_Volume; }
00256
00258 void SetFunctionId( const unsigned int& iFid )
00259 { this->m_FunctionId = iFid; }
00260
00261 virtual void ReleaseGlobalDataPointer(void *GlobalData) const
00262 { delete (GlobalDataStruct *) GlobalData; }
00263
00264 virtual ScalarValueType ComputeCurvatureTerm(const NeighborhoodType &,
00265 const FloatOffsetType &, GlobalDataStruct *gd );
00266
00269 virtual ScalarValueType LaplacianSmoothingSpeed(
00270 const NeighborhoodType &,
00271 const FloatOffsetType &, GlobalDataStruct * = 0) const
00272 { return NumericTraits<ScalarValueType>::One; }
00273
00276 virtual ScalarValueType CurvatureSpeed(const NeighborhoodType &,
00277 const FloatOffsetType &, GlobalDataStruct * = 0
00278 ) const
00279 { return NumericTraits<ScalarValueType>::One; }
00280
00281 protected:
00282
00283 RegionBasedLevelSetFunction();
00284 virtual ~RegionBasedLevelSetFunction() {}
00285
00287 InputImageConstPointer m_InitialImage;
00288
00290 FeatureImageConstPointer m_FeatureImage;
00291
00292 SharedDataPointer m_SharedData;
00293 HeavisideFunctionConstPointer m_DomainFunction;
00294
00296 ScalarValueType m_AreaWeight;
00297
00299 ScalarValueType m_Lambda1;
00300
00302 ScalarValueType m_Lambda2;
00303
00305 ScalarValueType m_OverlapPenaltyWeight;
00306
00308 ScalarValueType m_VolumeMatchingWeight;
00309
00311 ScalarValueType m_Volume;
00312
00314 ScalarValueType m_CurvatureWeight;
00315
00317 ScalarValueType m_LaplacianSmoothingWeight;
00318
00319 unsigned int m_FunctionId;
00320
00321 std::slice x_slice[itkGetStaticConstMacro(ImageDimension)];
00322 ::size_t m_Center;
00323 ::size_t m_xStride[itkGetStaticConstMacro(ImageDimension)];
00324
00325 static double m_WaveDT;
00326 static double m_DT;
00327
00328 void ComputeHImage();
00329
00332 ScalarValueType ComputeGlobalTerm(
00333 const ScalarValueType& imagePixel,
00334 const InputIndexType& inputIndex );
00335
00340 virtual ScalarValueType ComputeInternalTerm(const FeaturePixelType& iValue,
00341 const FeatureIndexType& iIdx ) = 0;
00342
00348 virtual ScalarValueType ComputeExternalTerm(const FeaturePixelType& iValue,
00349 const FeatureIndexType& iIdx ) = 0;
00350
00355 virtual ScalarValueType ComputeOverlapParameters( const FeatureIndexType& featIndex,
00356 ScalarValueType& pr ) = 0;
00357
00362 ScalarValueType ComputeVolumeRegularizationTerm( );
00363
00374 ScalarValueType ComputeLaplacianTerm( const NeighborhoodType &,
00375 const FloatOffsetType &, GlobalDataStruct *gd );
00376
00379 ScalarValueType ComputeLaplacian( GlobalDataStruct *gd );
00380
00382 void ComputeHessian( const NeighborhoodType &it,
00383 GlobalDataStruct *globalData );
00384
00386 virtual void ComputeParameters() = 0;
00387
00390 virtual void UpdateSharedDataParameters() = 0;
00391
00392 private:
00393 RegionBasedLevelSetFunction(const Self&);
00394 void operator=(const Self&);
00395 };
00396
00397 }
00398
00399 #ifndef ITK_MANUAL_INSTANTIATION
00400 #include "itkRegionBasedLevelSetFunction.txx"
00401 #endif
00402
00403 #endif
00404