Level Sets Framework Refactoring

Level Sets Framework Refactoring Arnaud Gelas, Kishore Mosaliganti, Nicolas Rannou, Lydie Souhait, Sean Megason Boston 02/03/2011

Outline Goals Principles Status Traits Fast Marching Image / Mesh Stopping Criterion Constrained Topology Shortest Path Computation Isotropic / Anisotropic Level Sets / GPU Plan Requirements

Goal: Generic Level Set Framework Multi-level set support simultaneous evolution of level sets Multi-channel support Vector or tensor data segmentation Level set representation Mesh-based (unstructured), image, or parametric Terms used in the PDE Add/delete terms in the update equation Topological constraints Stopping criterion RMS, Iterations, target points

STATUS Discrete Level Set

Traits TInputDomain TNode TOutputDomain TSuperclass

Traits - Base Class template< class TInputDomain, class TNode, class TOutputDomain, class TSuperclass > class LevelSetTraits { public: typedef [...] class NodePair : public std::pair< NodeType, OutputPixelType > [...] };

Traits - Image Specialization template<unsigned int VDimension,class TInputPixel,typename TOutputPixel >class ImageLevelSetTraits :public LevelSetTraits< Image< TInputPixel, VDimension >, Index< VDimension >, Image< TOutputPixel, VDimension >, ImageToImageFilter< Image< TInputPixel, VDimension >, Image< TOutputPixel, VDimension > > >

Traits - Mesh Specialization template< unsigned int VDimension,typename TInputPixel,class TInputMeshTraits,typename TOutputPixel,class TOutputMeshTraits >class MeshLevelSetTraits :public LevelSetTraits< Mesh< TInputPixel, VDimension, TInputMeshTraits >, typename TInputMeshTraits::PointIdentifier, Mesh< TOutputPixel, VDimension, TOutputMeshTraits >, MeshToMeshFilter< Mesh< TInputPixel, VDimension, TInputMeshTraits >, Mesh< TOutputPixel, VDimension, TOutputMeshTraits > >

Fast Marching Code available at the following address: https://github.com/arnaudgelas/itkFastMarching Numbers: 36 tests Tested Fedora 13, 14 (64 bits) Ubuntu 10.10 (64 bits) Mac OS-X 10.5, 10.6 Coverage: 80.49%

Stopping Criterion - Base Class class StoppingCriterionBase : public Object { public: virtualbool IsStatisfied() const = 0; virtualconst std::string GetDescription() const = 0; };

Stopping Criterion - Examples Threshold on the current value Equivalent to the current implementation ofitk::FastMarchingImageFilter<> Reached Target Nodes (One, Some, All), with possible overshoot offset Equivalent to the current implementation ofitk::FastMarchingUpWindGradientImageFilter<>

Constrained Topology [Tustisson'10 - Insight Journal 778] Escher's Ants as Metaphor: Topological Marching for the Well-Composed, Genus Zero Crowd

Minimal Path Extraction [Mueller'08 - Insight Journal 213] Fast Marching Minimal Path Extraction in ITK

Anisotropic • Several possible schemes • which one the best? • make it easy to implement any of these methods Isotropic / Anisotropic Isotropic Can be solved using current implementation

Requirements Possible performance improvement UpdateNeighbors() calls 2 * ImageDimension UpdateValue() Thread Pool ?

Fast Marching - Process ? Integration Process? Should we struggle for its integration (backward compatibility) ? Should we struggle a second time when integrating new level sets framework?

Fast Marching - Process ? update software guide? When ? How ? Any constraint?

FUTURE WORK

Git repository • Discrete Representations • DomainTraits • Iterators • Dense • Term container • Propagation • Advection • Curvature • Chan & Vese energy • Multithread • Reinitialization • Stopping Criterion Plan

Plan • Discrete Representations • Sparse • Constrained Topology • Multithread • Real time algorithm [Shi] • Parametric Representations • Splines • RBF

Discrete Level Sets - simplified view (a) while( ! m_StoppingCriterion->IsSatisfied() )(b) {(c) for each level set ls_i in the level set container(d) {(e) for each nodes n_j in the domain of ls_i(f) { (g) for each term t_k in the term container (h) { (i) Compute Term Value t_k( n_j, ls_i ) (j) Compute Term Contribution for time step computation (k) }(l) Evaluate the updated level set function ( delta( ls_i( n_j) ) )(m) } (n) } (o) Compute time step from CFL Condition (p) for each level set ls_i in the level set container (q) { (r) Update the level set function ls_i (s) Reinitialize to signed distance function (if requested by user)(t) }(u) }

GPU Involvement - 1 (a) while( !m_StoppingCriterion->IsSatisfied() )(b) {(c) for each level set ls_i in the level set container(d) {(e) for each nodes in the domain of ls_i(f) { (g) for each term t_k in the term container (h) { Compute Term Value t_k( n_j, ls_i ) GPU implementation during pixel updates at (i): Pixel neighborhood in image and level set is copied to GPU memory Terms are evaluated in the GPU function Each term will have a CPU and GPU implementation A term factory will call the GPU implementation Advantages: Minimal changes in the current proposed design Drawbacks: very bad according to performance

GPU Involvement -2 (a)while ( !m_StoppingCriterion->IsSatisfied() )(b) {(c) for each level set ls_i in the level set container(d) {(e) for each nodes in the domain of ls_i(f) { (g) for each term t_k in the term container (h) { Compute Term Value t_k( n_j, ls_i ) Entire while loop iteration (a) in GPU Everything is copied inside the GPU memory Advantages: Fastest solution in terms of performance Downside: memory limitation of the GPU (<2 Gb) Code duplication: CPU and GPU Note: Copy b/w memory 4Gb/s

GPU Involvement -3 In the last scenario, the code nesting is different: (a) while( !m_StoppingCriterion->IsSatisfied() )(b) {(c) for each level set ls_i in the level set container(d) {(e) for each term in the term container(f) {(g) for each nodes in the domain of ls_i (h) { (i) Evaluate the updated level set function [ ... ] In this one the GPU Implementation will occur for the most nested for loop (g) Keep copying the level set and image in each iteration in the GPU Second most optimal implementation for GPU No code duplication

Questions and Comments ?

Level Sets Framework Refactoring