~ubuntu-branches/ubuntu/precise/insighttoolkit/precise : revision 5

3

Program: Insight Segmentation & Registration Toolkit

4

Module: $RCSfile: itkParallelSparseFieldLevelSetImageFilter.h,v $

5

Language: C++

6

Date: $Date: 2008-03-03 13:58:45 $

7

Version: $Revision: 1.16 $

6

Date: $Date: 2008-10-17 16:30:48 $

7

Version: $Revision: 1.17 $

8

9

10

See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

14

PURPOSE. See the above copyright notices for more information.

15

16

=========================================================================*/

17

#ifndef __itkParallelSparseFieldLevelSetImageFilter_h_

18

#define __itkParallelSparseFieldLevelSetImageFilter_h_

17

#ifndef __itkParallelSparseFieldLevelSetImageFilter_h

18

#define __itkParallelSparseFieldLevelSetImageFilter_h

19

20

#include <vector>

21

#include "itkFiniteDifferenceImageFilter.h"

30

namespace itk

31

{

32

33

/**

33

/** \class ParallelSparseFieldLevelSetNode

34

* A data structure used in the ParallelSparsefieldlevelsetimagefilter to construct

35

* lists of indicies and other values.

36

*/

38

class ParallelSparseFieldLevelSetNode

39

{

40

public:

41

TNodeIndexType m_Index;

42

float m_Value;

41

TNodeIndexType m_Index;

42

float m_Value;

43

ParallelSparseFieldLevelSetNode *Next;

44

ParallelSparseFieldLevelSetNode *Previous;

45

};

74

class ParallelSparseFieldCityBlockNeighborList

75

{

76

public:

77

typedef TNeighborhoodType NeighborhoodType;

77

typedef TNeighborhoodType NeighborhoodType;

78

typedef typename NeighborhoodType::OffsetType OffsetType;

79

typedef typename NeighborhoodType::RadiusType RadiusType;

80

81

itkStaticConstMacro(Dimension, unsigned int, NeighborhoodType::Dimension);

82

83

const RadiusType &GetRadius() const

84

{

84

{

85

return m_Radius;

86

}

86

}

87

88

const unsigned int &GetArrayIndex(unsigned int i) const

89

{

89

{

90

return m_ArrayIndex[i];

91

}

91

}

92

93

const OffsetType &GetNeighborhoodOffset(unsigned int i) const

94

{

94

{

95

return m_NeighborhoodOffset[i];

96

}

96

}

97

98

const unsigned int &GetSize() const

99

{

99

{

100

return m_Size;

101

}

101

}

102

103

unsigned int GetStride(unsigned int i)

104

{

104

{

105

return m_StrideTable[i];

106

}

106

}

107

108

ParallelSparseFieldCityBlockNeighborList();

109

110

~ParallelSparseFieldCityBlockNeighborList()

111

{

111

{

112

m_ArrayIndex.clear();

113

m_NeighborhoodOffset.clear();

114

}

114

}

115

116

void Print(std::ostream &os) const;

117

118

private:

119

char pad1[128];

120

unsigned int m_Size;

121

RadiusType m_Radius;

119

char m_Pad1[128];

120

unsigned int m_Size;

121

RadiusType m_Radius;

122

std::vector<unsigned int> m_ArrayIndex;

123

std::vector<OffsetType> m_NeighborhoodOffset;

124

125

/** An internal table for keeping track of stride lengths in a neighborhood,

126

* i.e. the memory offsets between pixels along each dimensional axis. */

127

unsigned int m_StrideTable[Dimension];

128

char pad2[128];

128

char m_Pad2[128];

129

};

130

131

/**

249

{

250

public:

251

/** Standard class typedefs */

252

typedef ParallelSparseFieldLevelSetImageFilter Self;

252

typedef ParallelSparseFieldLevelSetImageFilter Self;

253

typedef FiniteDifferenceImageFilter<TInputImage, TOutputImage> Superclass;

254

typedef SmartPointer<Self> Pointer;

255

typedef SmartPointer<const Self> ConstPointer;

254

typedef SmartPointer<Self> Pointer;

255

typedef SmartPointer<const Self> ConstPointer;

256

257

/**Typedefs from the superclass */

258

typedef typename Superclass::TimeStepType TimeStepType;

258

typedef typename Superclass::TimeStepType TimeStepType;

259

typedef typename Superclass::FiniteDifferenceFunctionType FiniteDifferenceFunctionType;

260

typedef typename Superclass::RadiusType RadiusType;

261

typedef typename Superclass::NeighborhoodScalesType NeighborhoodScalesType;

260

typedef typename Superclass::RadiusType RadiusType;

261

typedef typename Superclass::NeighborhoodScalesType NeighborhoodScalesType;

262

263

/** Method for creation through the object factory. */

264

itkNewMacro(Self);

267

itkTypeMacro(ParallelSparseFieldLevelSetImageFilter, FiniteDifferenceImageFilter);

268

269

/** Information derived from the image types. */

270

typedef TInputImage InputImageType;

271

typedef TOutputImage OutputImageType;

270

typedef TInputImage InputImageType;

271

typedef TOutputImage OutputImageType;

272

typedef typename OutputImageType::IndexType IndexType;

273

274

itkStaticConstMacro(ImageDimension, unsigned int, TOutputImage::ImageDimension);

286

287

/** A list type used in the algorithm. */

288

typedef SparseFieldLayer<LayerNodeType> LayerType;

289

typedef typename LayerType::Pointer LayerPointerType;

289

typedef typename LayerType::Pointer LayerPointerType;

290

291

/** A type for a list of LayerPointerTypes */

292

typedef std::vector<LayerPointerType> LayerListType;

315

itkGetMacro(IsoSurfaceValue, ValueType);

316

317

LayerPointerType GetActiveListForIndex (const IndexType index)

318

{

319

// get the 'z' value for the index

320

const unsigned int indexZ= index[m_SplitAxis];

321

// get the thread in whose region the index lies

322

const unsigned int ThreadNum= this->GetThreadNumber (indexZ);

323

// get the active list for that thread

324

return m_Data[ThreadNum].m_Layers[0];

325

}

318

{

319

// get the 'z' value for the index

320

const unsigned int indexZ= index[m_SplitAxis];

321

// get the thread in whose region the index lies

322

const unsigned int ThreadNum= this->GetThreadNumber (indexZ);

323

// get the active list for that thread

324

return m_Data[ThreadNum].m_Layers[0];

325

}

326

327

#ifdef ITK_USE_CONCEPT_CHECKING

328

/** Begin concept checking */

388

389

/** The number of layers to use in the sparse field. Sparse field will

390

* consist of m_NumberOfLayers layers on both sides of a single active layer.

391

* This active layer is the interface of interest, i.e. the zero level set.*/

391

* This active layer is the interface of interest, i.e. the zero

392

* level set. */

392

393

StatusType m_NumberOfLayers;

393

394

395

/** An image of status values used internally by the algorithm. */

448

449

/** Adjusts the values associated with all the index layers of the sparse

449

450

* field by propagating out one layer at a time from the active set. This

450

451

* method also takes care of deleting nodes from the layers which have been

451

* marked in the status image as having been moved to other layers.*/

452

* marked in the status image as having been moved to other layers. */

452

453

void PropagateAllLayerValues();

453

454

455

/** Adjusts the values in a single layer "to" using values in a neighboring

469

470

471

/** Each thread allocates and initializes the data it will use by itself.

471

472

* This maintains the memory locality of data w.r.t. the thread in a shared

472

* memory environment.*/

473

* memory environment. */

473

474

void ThreadedAllocateData(unsigned int ThreadId);

474

475

void ThreadedInitializeData(unsigned int ThreadId, const ThreadRegionType & ThreadRegion);

475

476

496

497

498

/** Structure for managing thread-specific data */

498

499

struct ParallelSparseFieldLevelSetThreadStruct

499

{

500

{

500

501

ParallelSparseFieldLevelSetImageFilter* Filter;

501

502

TimeStepType* TimeStepList;

502

503

bool* ValidTimeStepList;

503

504

TimeStepType TimeStep;

504

};

505

};

505

506

507

/** This method calculates the change and does the update, i.e. one iteration

507

508

* of this iterative solver. A barrier class is used to synchronize

513

514

/** This method allows a subclass to override the way in which updates to

514

515

* output values are applied during each iteration. The default simply

515

516

* follows the standard finite difference scheme of scaling the change by the

516

* timestep and adding to the value of the previous iteration.*/

517

* timestep and adding to the value of the previous iteration. */

517

518

inline virtual ValueType ThreadedCalculateUpdateValue(const unsigned int itkNotUsed(ThreadId),

518

519

const IndexType itkNotUsed(index),

519

520

const TimeStepType &dt,

520

521

const ValueType &value,

521

522

const ValueType &change)

522

{

523

{

523

524

return (value + dt * change);

524

}

525

}

525

526

527

// This method can be overridden in derived classes.

527

528

// The pixel at 'index' is entering the active layer for thread 'ThreadId'.

528

// The outputimage at 'index' will have the value as given by the 'value' parameter.

529

// The outputimage at 'index' will have the value as given by the

530

// 'value' parameter.

529

531

virtual void ThreadedProcessPixelEnteringActiveLayer (const IndexType itkNotUsed(index),

530

532

const ValueType itkNotUsed(value),

531

533

const unsigned int itkNotUsed(ThreadId));

539

541

540

542

/** This method is not implemented or necessary for this solver */

541

543

TimeStepType CalculateChange()

542

{

544

{

543

545

return NumericTraits<TimeStepType>::Zero;

544

}

546

}

545

547

546

548

/** This method does the actual work of calculating change over a region

547

549

* supplied by the multithreading mechanism. */

680

682

681

683

/** Local data for each individual thread. */

682

684

struct ThreadData

683

{

685

{

684

686

char pad1 [128];

685

687

686

688

TimeStepType TimeStep;

719

721

/** Indicates whether to use m_Semaphore[0] or m_Semaphore[1] for signalling/waiting */

720

722

unsigned int m_SemaphoreArrayNumber;

721

723

722

char pad2 [128];

723

};

724

char m_Pad2 [128];

725

};

724

726

725

727

/** An array storing the individual (local) data structures for each thread. */

726

728

ThreadData *m_Data;