00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef __itkMultiphaseFiniteDifferenceImageFilter_h
00019 #define __itkMultiphaseFiniteDifferenceImageFilter_h
00020
00021 #include "itkArray.h"
00022 #include "itkInPlaceImageFilter.h"
00023 #include "itkFiniteDifferenceFunction.h"
00024 #include <vnl/vnl_vector.h>
00025 #include "itkImageRegionIterator.h"
00026
00027 #include "itkVector.h"
00028 #include "itkListSample.h"
00029 #include "itkKdTree.h"
00030 #include "itkKdTreeGenerator.h"
00031
00032 namespace itk {
00033
00161 template < class TInputImage,
00162 class TOutputImage,
00163 class TFiniteDifferenceFunction = FiniteDifferenceFunction<TOutputImage>,
00164 typename TIdCell = unsigned int >
00165 class ITK_EXPORT MultiphaseFiniteDifferenceImageFilter
00166 : public InPlaceImageFilter< TInputImage, TOutputImage >
00167 {
00168 public:
00169
00171 typedef MultiphaseFiniteDifferenceImageFilter Self;
00172 typedef InPlaceImageFilter<TInputImage, TOutputImage> Superclass;
00173 typedef SmartPointer<Self> Pointer;
00174 typedef SmartPointer<const Self> ConstPointer;
00175
00177 itkTypeMacro( MultiphaseFiniteDifferenceImageFilter, InPlaceImageFilter );
00178
00180 typedef TInputImage InputImageType;
00181 typedef typename InputImageType::Pointer InputImagePointer;
00182 typedef typename InputImageType::RegionType InputRegionType;
00183 typedef typename InputImageType::SizeType InputSizeType;
00184 typedef typename InputImageType::SpacingType InputSpacingType;
00185 typedef typename InputImageType::PointType InputPointType;
00186 typedef typename InputImageType::PixelType InputPixelType;
00187
00188 typedef TOutputImage OutputImageType;
00189 typedef typename OutputImageType::Pointer OutputImagePointer;
00190 typedef typename OutputImageType::PixelType OutputPixelType;
00191 typedef typename OutputImageType::RegionType OutputRegionType;
00192 typedef typename OutputImageType::SizeType OutputSizeType;
00193 typedef typename OutputImageType::SizeValueType OutputSizeValueType;
00194 typedef typename OutputImageType::IndexType OutputIndexType;
00195 typedef typename OutputImageType::IndexValueType OutputIndexValueType;
00196
00197 typedef TIdCell IdCellType;
00198 typedef std::vector< IdCellType > VectorIdCellType;
00199
00201 itkStaticConstMacro(ImageDimension, unsigned int, OutputImageType::ImageDimension);
00202
00206 typedef TFiniteDifferenceFunction FiniteDifferenceFunctionType;
00207 typedef typename FiniteDifferenceFunctionType::Pointer FiniteDifferenceFunctionPointer;
00208 typedef typename FiniteDifferenceFunctionType::TimeStepType TimeStepType;
00209 typedef typename std::vector< TimeStepType > TimeStepVectorType;
00210 typedef typename FiniteDifferenceFunctionType::RadiusType RadiusType;
00211
00212 typedef Vector< float, itkGetStaticConstMacro(ImageDimension) >
00213 CentroidVectorType;
00214 typedef Statistics::ListSample< CentroidVectorType > SampleType;
00215 typedef Statistics::KdTreeGenerator< SampleType > KdTreeGeneratorType;
00216 typedef typename KdTreeGeneratorType::Pointer KdTreeGeneratorPointer;
00217 typedef typename KdTreeGeneratorType::KdTreeType KdTreeType;
00218 typedef typename KdTreeType::Pointer KdTreePointer;
00219
00224 virtual const FiniteDifferenceFunctionPointer GetDifferenceFunction(
00225 const IdCellType& functionIndex ) const
00226 {
00227 if( functionIndex < m_FunctionCount )
00228 {
00229 return ( this->m_DifferenceFunctions[functionIndex] );
00230 }
00231 else
00232 {
00233 return 0;
00234 }
00235 }
00237
00242 virtual void SetDifferenceFunction( const IdCellType& functionIndex,
00243 FiniteDifferenceFunctionPointer function)
00244 {
00245 if( functionIndex < m_FunctionCount )
00246 {
00247 this->m_DifferenceFunctions[functionIndex] = function;
00248 }
00249 }
00251
00253 itkSetMacro(NumberOfIterations, unsigned int);
00254 itkGetConstReferenceMacro(NumberOfIterations, unsigned int);
00256
00259 itkSetMacro(UseImageSpacing,bool);
00260 itkBooleanMacro(UseImageSpacing);
00261 itkGetConstReferenceMacro(UseImageSpacing, bool);
00263
00266 itkSetMacro(MaximumRMSError, double);
00267 itkGetConstReferenceMacro(MaximumRMSError, double);
00269
00272 itkSetMacro(RMSChange, double);
00273 itkGetConstReferenceMacro(RMSChange, double);
00275
00277 itkSetMacro( InitializedState, bool );
00278 itkGetConstReferenceMacro( InitializedState, bool );
00279 itkBooleanMacro( InitializedState );
00281
00284 itkSetMacro( ManualReinitialization, bool );
00285 itkGetConstReferenceMacro( ManualReinitialization, bool );
00286 itkBooleanMacro( ManualReinitialization );
00288
00290 itkSetMacro( ElapsedIterations, unsigned int );
00291
00293 itkGetConstReferenceMacro( ElapsedIterations, unsigned int );
00294
00295 void SetLevelSet( const IdCellType & i, const InputImageType * levelSet )
00296 {
00297 m_LevelSet[i] = InputImageType::New();
00298 m_LevelSet[i]->SetRequestedRegion( levelSet->GetRequestedRegion() );
00299 m_LevelSet[i]->SetBufferedRegion( levelSet->GetBufferedRegion() );
00300 m_LevelSet[i]->SetLargestPossibleRegion( levelSet->GetLargestPossibleRegion() );
00301 m_LevelSet[i]->Allocate();
00302 m_LevelSet[i]->CopyInformation( levelSet );
00303
00304 ImageRegionConstIterator< InputImageType > in ( levelSet, levelSet->GetBufferedRegion());
00305 ImageRegionIterator< InputImageType > cp ( m_LevelSet[i], levelSet->GetBufferedRegion() );
00306
00307 in.GoToBegin();
00308 cp.GoToBegin();
00309
00310 while( !in.IsAtEnd() )
00311 {
00312 cp.Set( in.Get() );
00313 ++in;
00314 ++cp;
00315 }
00316 }
00317
00318 InputImagePointer GetLevelSet( const IdCellType& i )
00319 {
00320 if( i >= m_FunctionCount )
00321 {
00322 itkExceptionMacro("Request for level set #" << i
00323 << " but there are only " << m_FunctionCount );
00324 }
00325 else
00326 {
00327 return m_LevelSet[i];
00328 }
00329 }
00330
00331 void SetLookup ( VectorIdCellType lookup )
00332 {
00333 this->m_Lookup = lookup;
00334 }
00335
00336 void SetKdTree( KdTreeType * kdtree )
00337 {
00338 this->m_KdTree = kdtree;
00339 }
00340
00341 void SetFunctionCount( const IdCellType& n )
00342 {
00343 m_FunctionCount = n;
00344
00345 m_DifferenceFunctions.resize( m_FunctionCount, 0 );
00346
00347 RadiusType radius;
00348 radius.Fill( 1 );
00349
00350 for( unsigned int i = 0; i < this->m_FunctionCount; i++ )
00351 {
00352 this->m_DifferenceFunctions[i] = FiniteDifferenceFunctionType::New();
00353 this->m_DifferenceFunctions[i]->Initialize(radius);
00354 }
00355
00356
00357 m_LevelSet.resize( m_FunctionCount, 0 );
00358
00359
00360 this->m_Lookup.resize( m_FunctionCount );
00361
00362 IdCellType k = 1;
00363
00364 typedef typename std::vector< IdCellType >::iterator VectorIteratorType;
00365
00366 VectorIteratorType it = this->m_Lookup.begin();
00367
00368 while( it != this->m_Lookup.end() )
00369 {
00370 *it = k;
00371 ++it;
00372 ++k;
00373 }
00374 }
00375
00376 protected:
00377 MultiphaseFiniteDifferenceImageFilter()
00378 {
00379 this->m_KdTree = 0;
00380 this->m_ElapsedIterations = 0;
00381 this->m_MaximumRMSError = vnl_math::eps;
00382 this->m_RMSChange = 100.0;
00383 this->m_UseImageSpacing = true;
00384 this->m_ManualReinitialization = false;
00385 this->m_InitializedState = false;
00386 this->m_NumberOfIterations = NumericTraits<unsigned int>::max();
00387 this->m_FunctionCount = 0;
00388 this->InPlaceOn();
00389 }
00390
00391 ~MultiphaseFiniteDifferenceImageFilter(){}
00392
00393 IdCellType m_FunctionCount;
00394 std::vector< InputImagePointer > m_LevelSet;
00395 VectorIdCellType m_Lookup;
00396 KdTreePointer m_KdTree;
00397
00398 unsigned int m_ElapsedIterations;
00399 double m_MaximumRMSError;
00400 double m_RMSChange;
00401 unsigned int m_NumberOfIterations;
00402
00404 std::vector< FiniteDifferenceFunctionPointer > m_DifferenceFunctions;
00405
00408 bool m_UseImageSpacing;
00409
00410 void PrintSelf(std::ostream& os, Indent indent) const;
00411
00413 virtual void AllocateUpdateBuffer() = 0;
00414
00418 virtual void ApplyUpdate(TimeStepType dt) = 0;
00419
00425 virtual TimeStepType CalculateChange() = 0;
00426
00430 virtual void CopyInputToOutput() = 0;
00431
00435 virtual void GenerateData();
00436
00448 virtual void GenerateInputRequestedRegion();
00449
00452 virtual bool Halt();
00453
00463 virtual bool ThreadedHalt(void *itkNotUsed(threadInfo))
00464 {
00465 return this->Halt();
00466 }
00467
00473 virtual void Initialize() { };
00474
00481 virtual void InitializeIteration()
00482 {
00483 for( IdCellType i = 0; i < this->m_FunctionCount; i++ )
00484 {
00485 this->m_DifferenceFunctions[i]->InitializeIteration();
00486 }
00487 }
00489
00503 inline TimeStepType ResolveTimeStep(const TimeStepVectorType& timeStepList,
00504 const std::vector< bool >& valid );
00505
00508 virtual void PostProcessOutput() {}
00509
00510
00511 private:
00512 MultiphaseFiniteDifferenceImageFilter(const Self&);
00513
00514 void operator=(const Self&);
00515
00518 bool m_ManualReinitialization;
00519
00521 bool m_InitializedState;
00522 };
00523
00524 }
00525
00526 #ifndef ITK_MANUAL_INSTANTIATION
00527 #include "itkMultiphaseFiniteDifferenceImageFilter.txx"
00528 #endif
00529
00530 #endif
00531