diff --git a/gems/CMakeLists.txt b/gems/CMakeLists.txt
index 3a62a57..c6f5db1 100755
--- a/gems/CMakeLists.txt
+++ b/gems/CMakeLists.txt
@@ -160,13 +160,13 @@ set(SOURCES_DYN
 
 # gems libary
 add_library(kvlGEMSCommon ${SOURCES})
-target_link_libraries(kvlGEMSCommon ${ITK_LIBRARIES} ${ZLIB_LIBRARIES})
+target_link_libraries(kvlGEMSCommon ${ZLIB_LIBRARIES} ${ITK_LIBRARIES})
 
 # gems library with dynamic meshes.  Needed for `kvlBuildAtlasMesh` but causes
 # samseg to run slower, so only link against `kvlGEMSCommon_dynmesh` when nessesary.
 add_library(kvlGEMSCommon_dynmesh ${SOURCES_DYN})
 target_compile_definitions(kvlGEMSCommon_dynmesh PRIVATE -DUSE_DYNAMIC_MESH)
-target_link_libraries(kvlGEMSCommon_dynmesh ${ITK_LIBRARIES} ${ZLIB_LIBRARIES})
+target_link_libraries(kvlGEMSCommon_dynmesh ${ZLIB_LIBRARIES} ${ITK_LIBRARIES})
 
 ## !!!BAD!!! both `kvlGEMSCommon_dynmesh` `kvlGEMSCommon` will have the same md5sum
 ## Useful for testing though
diff --git a/gems/Executables/kvlAtlasMeshBuilder.cxx b/gems/Executables/kvlAtlasMeshBuilder.cxx
index ff113ee..98f279b 100755
--- a/gems/Executables/kvlAtlasMeshBuilder.cxx
+++ b/gems/Executables/kvlAtlasMeshBuilder.cxx
@@ -282,7 +282,7 @@ ::TryToCollapseFast( const AtlasMeshCollection* miniCollection,
   AtlasMeshCollection::Pointer  child;
   if ( !miniCollection->GetCollapsed( edgeId, child, disappearingCells, unifiedVertexId ) )
     {
-    return 0;
+    return nullptr;
     }
 
   const AtlasMesh::PointIdentifier  unifiedPointId =
@@ -405,7 +405,7 @@ ::TryToCollapseFast( const AtlasMeshCollection* miniCollection,
         }
       else
         {
-        return 0;
+        return nullptr;
         }
 
       }  // End test if point can move
@@ -418,7 +418,7 @@ ::TryToCollapseFast( const AtlasMeshCollection* miniCollection,
   if ( !costFunction->GetValue( optimalPosition, miniDataCost, miniAlphasCost, miniPositionCost ) )
     {
     //std::cout << "not possible" << std::endl;
-    return 0;
+    return nullptr;
     }
   //std::cout << "       miniDataCost     : " << miniDataCost << std::endl;
   //std::cout << "       miniAlphasCost   : " << miniAlphasCost << std::endl;
@@ -429,7 +429,7 @@ ::TryToCollapseFast( const AtlasMeshCollection* miniCollection,
   if ( std::isnan( miniDataCost + miniAlphasCost + miniPositionCost ) ||
        std::isinf( miniDataCost + miniAlphasCost + miniPositionCost ) )
     {
-    return 0;
+    return nullptr;
     }
   else
     {
@@ -496,11 +496,11 @@ ::TryToRetainFast( const AtlasMeshCollection* miniCollectionConst,
   // TODO: should be iterated!
   if ( !this->OptimizeReferencePositionFast( miniCollection, vertex0Id, false ) )
     {
-    return 0;
+    return nullptr;
     }
   if ( !this->OptimizeReferencePositionFast( miniCollection, vertex1Id, false ) )
     {
-    return 0;
+    return nullptr;
     }
 
 
@@ -513,7 +513,7 @@ ::TryToRetainFast( const AtlasMeshCollection* miniCollectionConst,
   if ( ( std::isnan( miniDataCost + miniAlphasCost + miniPositionCost ) ) ||
        ( std::isinf( miniDataCost + miniAlphasCost + miniPositionCost ) ) )
     {
-    return 0;
+    return nullptr;
     }
   else
     {
@@ -1304,16 +1304,29 @@ ::GetRandomizedEdgesAsSet() const
 //
 //
 //
+#if ITK_VERSION_MAJOR >= 5
+itk::ITK_THREAD_RETURN_TYPE
+AtlasMeshBuilder
+::LoadBalancedThreaderCallback( void *arg )
+#else  
 ITK_THREAD_RETURN_TYPE
 AtlasMeshBuilder
 ::LoadBalancedThreaderCallback( void *arg )
+#endif    
 {
 
   // Retrieve the input arguments
+#if ITK_VERSION_MAJOR >= 5
+  const int  threadId = ((itk::MultiThreaderBase::WorkUnitInfo *)(arg))->WorkUnitID;
+  //const int  threadCount = ((itk::MultiThreaderBase::WorkUnitInfo *)(arg))->NumberOfWorkUnits;
+
+  LoadBalancedThreadStruct*  str = (LoadBalancedThreadStruct *)(((itk::MultiThreaderBase::WorkUnitInfo *)(arg))->UserData);  
+#else  
   const int  threadId = ((itk::MultiThreader::ThreadInfoStruct *)(arg))->ThreadID;
   //const int  threadCount = ((itk::MultiThreader::ThreadInfoStruct *)(arg))->NumberOfThreads;
 
   LoadBalancedThreadStruct*  str = (LoadBalancedThreadStruct *)(((itk::MultiThreader::ThreadInfoStruct *)(arg))->UserData);
+#endif  
 
   const int  numberOfEdgesToAnalyze = str->m_Edges.size();
   int  numberOfEdgesAnalyzed = 0;
@@ -1358,7 +1371,11 @@ ::LoadBalancedThreaderCallback( void *arg )
       }  
     }
 
+#if ITK_VERSION_MAJOR >= 5
+  return itk::ITK_THREAD_RETURN_DEFAULT_VALUE;
+#else    
   return ITK_THREAD_RETURN_VALUE;
+#endif  
 }
 
 
diff --git a/gems/Executables/kvlAtlasMeshBuilder.h b/gems/Executables/kvlAtlasMeshBuilder.h
index d044dc2..52baea1 100755
--- a/gems/Executables/kvlAtlasMeshBuilder.h
+++ b/gems/Executables/kvlAtlasMeshBuilder.h
@@ -1,6 +1,10 @@
 #ifndef __kvlAtlasMeshBuilder_h
 #define __kvlAtlasMeshBuilder_h
 
+#if ITK_VERSION_MAJOR >= 5
+#include <mutex>
+#endif
+
 #include "kvlMultiResolutionAtlasMesher.h"
 #include "vnl/vnl_sample.h"
 #include "itkTimeProbe.h"
@@ -11,17 +15,29 @@ namespace kvl
 {
 
 
+#if ITK_VERSION_MAJOR >= 5
+class AtlasMeshBuilderMutexLock: public std::mutex
+#else
 class AtlasMeshBuilderMutexLock: public itk::SimpleFastMutexLock
+#endif
 {
 public:
   /** Standard class typedefs.  */
   typedef AtlasMeshBuilderMutexLock       Self;
+#if ITK_VERSION_MAJOR >= 5
+  typedef std::mutex                      Superclass;
+#else  
   typedef itk::SimpleFastMutexLock        Superclass;
+#endif  
 
   /** Lock access. */
   void DescriptiveLock( const std::string& description )
     {
+#if ITK_VERSION_MAJOR >= 5
+    Superclass::lock();      
+#else      
     Superclass::Lock();
+#endif    
     m_TimeProbe.Start();
     m_Description = description;
     }
@@ -31,7 +47,11 @@ class AtlasMeshBuilderMutexLock: public itk::SimpleFastMutexLock
     {
     m_TimeProbe.Stop();
     std::cout << m_Description << ": unlocking mutex after " << m_TimeProbe.GetMean() << " seconds" << std::endl;
+#if ITK_VERSION_MAJOR >= 5
+    Superclass::unlock();
+#else
     Superclass::Unlock();
+#endif    
     }
 
 protected:
@@ -58,7 +78,11 @@ class AtlasMeshBuilderHelper
   AtlasMeshBuilderHelper( AtlasMeshBuilderMutexLock&  mutex, std::map< AtlasMesh::PointIdentifier, int >&  pointOccupancies )
   : m_Mutex( mutex ), m_MutexIsLocked( true ), m_PointOccupancies( pointOccupancies )
     {
-    m_Mutex.Lock();
+#if ITK_VERSION_MAJOR >= 5
+      m_Mutex.lock();
+#else      
+      m_Mutex.Lock();
+#endif
     }
 
 
@@ -74,7 +98,11 @@ class AtlasMeshBuilderHelper
     if ( !m_MutexIsLocked )
       {
       m_MutexIsLocked = true;
+#if ITK_VERSION_MAJOR >= 5
+      m_Mutex.lock();
+#else      
       m_Mutex.Lock();
+#endif      
       }
       
     }
@@ -84,7 +112,11 @@ class AtlasMeshBuilderHelper
     if ( m_MutexIsLocked )
       {
       m_MutexIsLocked = false;
+#if ITK_VERSION_MAJOR >= 5
+      m_Mutex.unlock();
+#else      
       m_Mutex.Unlock();
+#endif      
       }
     }  
     
@@ -319,7 +351,11 @@ protected :
   /** Static function used as a "callback" by the MultiThreader.  The threading
    * library will call this routine for each thread, which will delegate the
    * control to ThreadedGenerateData(). */
+#if ITK_VERSION_MAJOR >= 5
+  static itk::ITK_THREAD_RETURN_TYPE LoadBalancedThreaderCallback( void *arg );
+#else  
   static ITK_THREAD_RETURN_TYPE LoadBalancedThreaderCallback( void *arg );
+#endif
 
   /** Internal structure used for passing image data into the threading library */
   struct LoadBalancedThreadStruct
diff --git a/gems/Executables/kvlBuildAtlasMesh.cxx b/gems/Executables/kvlBuildAtlasMesh.cxx
index 00aecdd..ebf55b0 100755
--- a/gems/Executables/kvlBuildAtlasMesh.cxx
+++ b/gems/Executables/kvlBuildAtlasMesh.cxx
@@ -262,7 +262,7 @@ int main( int argc, char** argv )
 
 
   // If explicitStartCollection exists in the current directory, use it
-  kvl::AtlasMeshCollection::Pointer  explicitStartCollection = 0;
+  kvl::AtlasMeshCollection::Pointer  explicitStartCollection = nullptr;
   const std::string  explicitStartCollectionFileName = "explicitStartCollection.gz";
   if ( itksys::SystemTools::FileExists( explicitStartCollectionFileName.c_str(), true ) )
     {
diff --git a/gems/itkMGHImageIO.h b/gems/itkMGHImageIO.h
index 2bf9b88..ad92ebe 100755
--- a/gems/itkMGHImageIO.h
+++ b/gems/itkMGHImageIO.h
@@ -8,7 +8,11 @@
 // ITK includes
 #include "itkImageIOBase.h"
 #include "itkIOCommon.h"
+#if ITK_VERSION_MAJOR >= 5
+#include "itkMacro.h"
+#else
 #include "itkExceptionObject.h"
+#endif
 #include "itkByteSwapper.h"
 #include "itkMetaDataObject.h"
 #include "itkMatrix.h"
diff --git a/gems/kvlAtlasMeshCollection.cxx b/gems/kvlAtlasMeshCollection.cxx
index 8c5aff6..76a63f4 100755
--- a/gems/kvlAtlasMeshCollection.cxx
+++ b/gems/kvlAtlasMeshCollection.cxx
@@ -1590,7 +1590,7 @@ ::GetCollapsed( AtlasMesh::CellIdentifier edgeId,
   // disappears, and first point, which gets the previously determined new position
   //std::cout << "Creating positions" << std::endl;
   std::vector< PointsContainerType::Pointer >  collapsedPositions;
-  PointsContainerType::Pointer   collapsedReferencePosition = 0;
+  PointsContainerType::Pointer   collapsedReferencePosition = nullptr;
   for ( unsigned int meshNumber = 0; meshNumber < this->GetNumberOfMeshes()+1; meshNumber++ )
     {
     PointsContainerType::Pointer  thisPosition;
@@ -1857,7 +1857,7 @@ ::GetRegionGrown( AtlasMesh::CellIdentifier seedId, unsigned int radius, bool ma
   if ( !( m_Cells->IndexExists( seedId ) ) )
     {
     // Cell simply doesn't exist.
-    return 0;
+    return nullptr;
     }
  
   // Get the cell
@@ -2353,7 +2353,7 @@ ::GetUpsampled() const
 
   // Loop over all meshes, including reference
   std::vector< PointsContainerType::Pointer >  upsampledPositions;
-  PointsContainerType::Pointer   upsampledReferencePosition = 0;
+  PointsContainerType::Pointer   upsampledReferencePosition = nullptr;
   CellsContainerType::Pointer  upsampledCells;
   for ( unsigned int meshNumber=0; meshNumber < this->GetNumberOfMeshes()+1; meshNumber++ )
     {
@@ -2985,7 +2985,7 @@ ::GetUpsampled() const
     if ( upsampledPositions[ meshNumber ]->Size() != upsampledReferencePosition->Size() )
       {
       std::cerr << "Upsampling failed because of numerical inaccuracies!" << std::endl;
-      return 0;
+      return nullptr;
       }
     }
       
@@ -3085,7 +3085,7 @@ ::GetEdgeSplitted( AtlasMesh::CellIdentifier edgeId,
 
   return this->GetEdgeSplitted( edgeId, newVertexId, newPointId, transverseEdge0IdDummy, transverseEdge1IdDummy );
 #else
-  return 0;
+  return nullptr;
 #endif
 
 }
@@ -3112,7 +3112,7 @@ ::GetEdgeSwapped( AtlasMesh::CellIdentifier edgeId ) const
   AtlasMesh::CellIdentifier  newEdgeIdDummy;
   return this->GetEdgeSwapped( edgeId, newVertexId, newPointId, newEdgeIdDummy );
 #else
-  return 0;
+  return nullptr;
 #endif
 }
 
@@ -3232,7 +3232,7 @@ ::GetEdgeSwapped( AtlasMesh::CellIdentifier edgeId,
   return collapsed;
 
 #else
-  return 0;
+  return nullptr;
 
 #endif
 }                 
@@ -3591,7 +3591,7 @@ ::GetEdgeSplitted( AtlasMesh::CellIdentifier  edgeId,
   return splitted;
 
 #else
-  return 0;
+  return nullptr;
 #endif
 
 }
diff --git a/gems/kvlAtlasMeshDeformationConjugateGradientOptimizer.cxx b/gems/kvlAtlasMeshDeformationConjugateGradientOptimizer.cxx
index 388dd46..bb61a6d 100755
--- a/gems/kvlAtlasMeshDeformationConjugateGradientOptimizer.cxx
+++ b/gems/kvlAtlasMeshDeformationConjugateGradientOptimizer.cxx
@@ -60,7 +60,7 @@ ::FindAndOptimizeNewSearchDirection()
   // 
   // Part I: Decide on a new search direction
   //
-  AtlasPositionGradientContainerType::Pointer  searchDirection = 0;
+  AtlasPositionGradientContainerType::Pointer  searchDirection = nullptr;
   bool  startingOrRestarting = false;
   if ( this->GetIterationNumber() == 0 )
     {
diff --git a/gems/kvlAtlasMeshDeformationFixedStepGradientDescentOptimizer.cxx b/gems/kvlAtlasMeshDeformationFixedStepGradientDescentOptimizer.cxx
index 8e46873..c9e5ab7 100755
--- a/gems/kvlAtlasMeshDeformationFixedStepGradientDescentOptimizer.cxx
+++ b/gems/kvlAtlasMeshDeformationFixedStepGradientDescentOptimizer.cxx
@@ -46,7 +46,7 @@ ::FindAndOptimizeNewSearchDirection()
   // Try to add the scaled gradient to the current position to obtain the trial position
   const double  alpha = -( m_StepSize / maximumGradientMagnitude );
   double  maximalDeformation = 0.0;
-  AtlasMesh::PointsContainer::Pointer  trialPosition = 0;
+  AtlasMesh::PointsContainer::Pointer  trialPosition = nullptr;
   this->AddDeformation( m_Position, alpha, m_Gradient, trialPosition, maximalDeformation );
   if ( m_Verbose )
     {
@@ -54,7 +54,7 @@ ::FindAndOptimizeNewSearchDirection()
     }
   
   // Try out this new position
-  AtlasPositionGradientContainerType::Pointer  trialGradient = 0;
+  AtlasPositionGradientContainerType::Pointer  trialGradient = nullptr;
   double  trialCost = 0.0;
   this->GetCostAndGradient( trialPosition, trialCost, trialGradient );
   if ( m_Verbose )
diff --git a/gems/kvlAtlasMeshDeformationOptimizer.cxx b/gems/kvlAtlasMeshDeformationOptimizer.cxx
index 6428ec4..2ed696e 100755
--- a/gems/kvlAtlasMeshDeformationOptimizer.cxx
+++ b/gems/kvlAtlasMeshDeformationOptimizer.cxx
@@ -376,11 +376,11 @@ ::DoLineSearch( const AtlasMesh::PointsContainer*  startPosition,
   const int  maximumOfBracketingIterations = 10;
   double  lowAlpha = 0.0;
   double  lowCost = 0.0;
-  AtlasPositionGradientContainerType::ConstPointer  lowGradient = 0;
+  AtlasPositionGradientContainerType::ConstPointer  lowGradient = nullptr;
   double  lowDirectionalDerivative = 0.0;
   double  highAlpha = 0.0;
   double  highCost = 0.0;
-  AtlasPositionGradientContainerType::ConstPointer  highGradient = 0;
+  AtlasPositionGradientContainerType::ConstPointer  highGradient = nullptr;
   double  highDirectionalDerivative = 0.0;
       
   for ( int bracketingIterationNumber = 0; 
@@ -388,12 +388,12 @@ ::DoLineSearch( const AtlasMesh::PointsContainer*  startPosition,
         bracketingIterationNumber++ )
     {     
     // Evaluate current alpha: [ cost gradient ] = tryFunction( x + alpha * p );
-    AtlasMesh::PointsContainer::Pointer  position = 0;
+    AtlasMesh::PointsContainer::Pointer  position = nullptr;
     double  maximalDeformation = 0.0;
     this->AddDeformation( startPosition, alpha, searchDirection, 
                           position, maximalDeformation );
     double  cost;
-    AtlasPositionGradientContainerType::Pointer  gradient = 0;
+    AtlasPositionGradientContainerType::Pointer  gradient = nullptr;
     this->GetCostAndGradient( position, cost, gradient );
     const double  directionalDerivative = this->ComputeInnerProduct( gradient, searchDirection ); // gradient' * p
 
@@ -721,12 +721,12 @@ ::DoLineSearch( const AtlasMesh::PointsContainer*  startPosition,
 #endif
     
     // Evaluate cost function: [ cost gradient ] = tryFunction( x + alpha * p );
-    AtlasMesh::PointsContainer::Pointer  position = 0;
+    AtlasMesh::PointsContainer::Pointer  position = nullptr;
     double  maximalDeformation = 0.0;
     this->AddDeformation( startPosition, alpha, searchDirection, 
                           position, maximalDeformation );
     double  cost;
-    AtlasPositionGradientContainerType::Pointer  gradient = 0;
+    AtlasPositionGradientContainerType::Pointer  gradient = nullptr;
     this->GetCostAndGradient( position, cost, gradient );
     const double  directionalDerivative = this->ComputeInnerProduct( gradient, searchDirection ); // gradient' * p
 
diff --git a/gems/kvlAtlasMeshPositionCostAndGradientCalculator.cxx b/gems/kvlAtlasMeshPositionCostAndGradientCalculator.cxx
index caad092..7bcc299 100755
--- a/gems/kvlAtlasMeshPositionCostAndGradientCalculator.cxx
+++ b/gems/kvlAtlasMeshPositionCostAndGradientCalculator.cxx
@@ -250,7 +250,7 @@ AtlasMeshPositionCostAndGradientCalculator
 ::SetMeshToImageTransform( const TransformType* meshToImageTransform )
 {
   
-  TransformType::Pointer  inScopeHolder = 0;
+  TransformType::Pointer  inScopeHolder = nullptr;
   if ( !meshToImageTransform )
     {
     //meshToImageTransform = TransformType::New();
diff --git a/gems/kvlAtlasMeshRasterizor.cxx b/gems/kvlAtlasMeshRasterizor.cxx
index c7c1c56..e3d4e13 100755
--- a/gems/kvlAtlasMeshRasterizor.cxx
+++ b/gems/kvlAtlasMeshRasterizor.cxx
@@ -1,6 +1,12 @@
 #include "kvlAtlasMeshRasterizor.h"
 
+#if ITK_VERSION_MAJOR >= 5
+#include <itkMultiThreaderBase.h>
+#include <mutex>
+static std::mutex rasterizorMutex;
+#else
 static itk::SimpleFastMutexLock rasterizorMutex;
+#endif
 
 
 
@@ -14,7 +20,11 @@ namespace kvl
 AtlasMeshRasterizor
 ::AtlasMeshRasterizor()
 {
+#if ITK_VERSION_MAJOR >= 5
+  m_NumberOfThreads = itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads();
+#else  
   m_NumberOfThreads = itk::MultiThreader::GetGlobalDefaultNumberOfThreads();
+#endif  
 }
 
 
@@ -42,9 +52,14 @@ ::Rasterize( const AtlasMesh* mesh )
     }
 
   // Set up the multithreader
+#if ITK_VERSION_MAJOR >= 5
+  itk::MultiThreaderBase::Pointer  threader = itk::MultiThreaderBase::New();
+  threader->SetNumberOfWorkUnits( this->GetNumberOfThreads() );
+#else
   itk::MultiThreader::Pointer  threader = itk::MultiThreader::New();
   threader->SetNumberOfThreads( this->GetNumberOfThreads() );
   //threader->SetNumberOfThreads( 1 );
+#endif  
   threader->SetSingleMethod( this->ThreaderCallback, &str );
 
   // Let the beast go
@@ -59,15 +74,27 @@ ::Rasterize( const AtlasMesh* mesh )
 //
 //
 //
+#if ITK_VERSION_MAJOR >= 5
+itk::ITK_THREAD_RETURN_TYPE
+AtlasMeshRasterizor
+::ThreaderCallback( void *arg )
+#else  
 ITK_THREAD_RETURN_TYPE
 AtlasMeshRasterizor
 ::ThreaderCallback( void *arg )
+#endif  
 {
 
   // Retrieve the input arguments
+#if ITK_VERSION_MAJOR >= 5
+  const int  threadNumber = ((itk::MultiThreaderBase::WorkUnitInfo *)(arg))->WorkUnitID;
+  const int  numberOfThreads = ((itk::MultiThreaderBase::WorkUnitInfo *)(arg))->NumberOfWorkUnits;
+  ThreadStruct*  str = (ThreadStruct *)(((itk::MultiThreaderBase::WorkUnitInfo *)(arg))->UserData);
+#else  
   const int  threadNumber = ((itk::MultiThreader::ThreadInfoStruct *)(arg))->ThreadID;
   const int  numberOfThreads = ((itk::MultiThreader::ThreadInfoStruct *)(arg))->NumberOfThreads;
   ThreadStruct*  str = (ThreadStruct *)(((itk::MultiThreader::ThreadInfoStruct *)(arg))->UserData);
+#endif  
 
   
 #if 1  
@@ -116,7 +143,11 @@ ::ThreaderCallback( void *arg )
     if ( str->m_TetrahedronIds.size() == 0 )
       {
       rasterizorMutex.Unlock();
+#if ITK_VERSION_MAJOR >= 5
+      return itk::ITK_THREAD_RETURN_DEFAULT_VALUE;
+#else        
       return ITK_THREAD_RETURN_VALUE;
+#endif      
       }
 
     // Let's define how many tetrahedra this thread is going to take on 
@@ -151,7 +182,11 @@ ::ThreaderCallback( void *arg )
         str->m_TetrahedronIds.clear();
         rasterizorMutex.Unlock();
           
+#if ITK_VERSION_MAJOR >= 5
+        return itk::ITK_THREAD_RETURN_DEFAULT_VALUE;
+#else  	
         return ITK_THREAD_RETURN_VALUE;
+#endif	
         }
         
       }  
@@ -161,7 +196,11 @@ ::ThreaderCallback( void *arg )
 #endif
     
   
+#if ITK_VERSION_MAJOR >= 5
+  return itk::ITK_THREAD_RETURN_DEFAULT_VALUE;
+#else  
   return ITK_THREAD_RETURN_VALUE;
+#endif  
 }
 
 
diff --git a/gems/kvlAtlasMeshRasterizor.h b/gems/kvlAtlasMeshRasterizor.h
index 1d7f658..d412a98 100755
--- a/gems/kvlAtlasMeshRasterizor.h
+++ b/gems/kvlAtlasMeshRasterizor.h
@@ -68,7 +68,11 @@ class AtlasMeshRasterizor : public itk::Object
   /** Static function used as a "callback" by the MultiThreader.  The threading
    * library will call this routine for each thread, which will delegate the
    * control to ThreadedGenerateData(). */
+#if ITK_VERSION_MAJOR >= 5
+  static itk::ITK_THREAD_RETURN_TYPE ThreaderCallback( void *arg );
+#else  
   static ITK_THREAD_RETURN_TYPE ThreaderCallback( void *arg );
+#endif  
   
   /** Internal structure used for passing information to the threading library */
   struct ThreadStruct
diff --git a/gems/kvlAtlasMeshSmoother.cxx b/gems/kvlAtlasMeshSmoother.cxx
index e9f2a4e..d1c0214 100755
--- a/gems/kvlAtlasMeshSmoother.cxx
+++ b/gems/kvlAtlasMeshSmoother.cxx
@@ -64,7 +64,7 @@ ::GetSmoothedMeshCollection()
   //  std::cout << "Smoothing mesh collection with sigma = " << m_Sigma << "..." << std::flush;
 
   // Construct container to hold smoothed alphas
-  AtlasMesh::PointDataContainer::Pointer  smoothedParameters = 0;
+  AtlasMesh::PointDataContainer::Pointer  smoothedParameters = nullptr;
 
   
   if ( ( m_Sigma0 == 0 ) && ( m_Sigma1 == 0 ) && ( m_Sigma2 ==0 ) )
diff --git a/gems/kvlAtlasMeshValueDrawer.h b/gems/kvlAtlasMeshValueDrawer.h
index ac98db7..9b5df10 100755
--- a/gems/kvlAtlasMeshValueDrawer.h
+++ b/gems/kvlAtlasMeshValueDrawer.h
@@ -60,7 +60,7 @@ public :
   
   int m_NumFrames = 0;
   double const * m_Values = 0;
-  ImageType::Pointer m_Image = 0;
+  ImageType::Pointer m_Image = nullptr;
   
 };
 
diff --git a/gems/kvlAtlasParameterEstimator.cxx b/gems/kvlAtlasParameterEstimator.cxx
index 1b0b1cb..35de249 100755
--- a/gems/kvlAtlasParameterEstimator.cxx
+++ b/gems/kvlAtlasParameterEstimator.cxx
@@ -1,5 +1,9 @@
 #include "kvlAtlasParameterEstimator.h"
 
+#if ITK_VERSION_MAJOR >= 5
+#include <itkMultiThreaderBase.h>
+#endif
+
 #include "itkImageRegionConstIterator.h"
 #include "kvlAtlasMeshLabelImageStatisticsCollector.h"
 #include "kvlAtlasMeshToLabelImageCostAndGradientCalculator.h"
@@ -40,7 +44,11 @@ ::AtlasParameterEstimator()
 
   m_PositionOptimizer = LBFGS;
 
+#if ITK_VERSION_MAJOR >= 5
+  m_NumberOfThreads = itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads();
+#else  
   m_NumberOfThreads = itk::MultiThreader::GetGlobalDefaultNumberOfThreads();
+#endif  
   
 }
 
@@ -363,7 +371,7 @@ ::EstimatePosition( unsigned int  labelImageNumber )
   
   
   // Set up a deformation optimizer and pass the gradient calculator on to it
-  AtlasMeshDeformationOptimizer::Pointer  optimizer = 0;
+  AtlasMeshDeformationOptimizer::Pointer  optimizer = nullptr;
   switch( m_PositionOptimizer ) 
     {
     case FIXED_STEP_GRADIENT_DESCENT: 
@@ -430,7 +438,7 @@ ::CalculateCurrentPositionCostAndGradient( unsigned int labelImageNumber, double
   if ( ( labelImageNumber >= m_MeshCollection->GetNumberOfMeshes() ) || 
        ( labelImageNumber >= m_LabelImages.size() ) )
     {
-    return 0;
+    return nullptr;
     }    
 
   // Set up gradient calculator
diff --git a/gems/kvlDSWbetaMMLikelihoodImageFilter.hxx b/gems/kvlDSWbetaMMLikelihoodImageFilter.hxx
index 0c6e4f0..1f586a8 100644
--- a/gems/kvlDSWbetaMMLikelihoodImageFilter.hxx
+++ b/gems/kvlDSWbetaMMLikelihoodImageFilter.hxx
@@ -16,6 +16,11 @@ DSWbetaMMLikelihoodImageFilter< TInputImage >
 ::DSWbetaMMLikelihoodImageFilter()
 {
     //m_likelihoodFilter = GMMLikelihoodFilterType::New();
+#if ITK_VERSION_MAJOR >= 5 
+  // use classic void ThreadedGenerateData( const OutputRegionType& threadRegion, ThreadIdType threadId )
+  // instead of default new signature void DynamicThreadedGenerateData( const OutputRegionType& threadRegion )
+  this->DynamicMultiThreadingOff();
+#endif  
 }
 
 
diff --git a/gems/kvlFrobMMLikelihoodImageFilter.hxx b/gems/kvlFrobMMLikelihoodImageFilter.hxx
index 001ccb6..ea941a5 100644
--- a/gems/kvlFrobMMLikelihoodImageFilter.hxx
+++ b/gems/kvlFrobMMLikelihoodImageFilter.hxx
@@ -16,6 +16,11 @@ FrobMMLikelihoodImageFilter< TInputImage >
 ::FrobMMLikelihoodImageFilter()
 {
     //m_likelihoodFilter = GMMLikelihoodFilterType::New();
+#if ITK_VERSION_MAJOR >= 5 
+  // use classic void ThreadedGenerateData( const OutputRegionType& threadRegion, ThreadIdType threadId )
+  // instead of default new signature void DynamicThreadedGenerateData( const OutputRegionType& threadRegion )
+  this->DynamicMultiThreadingOff();
+#endif 
 }
 
 
diff --git a/gems/kvlGMMLikelihoodImageFilter.hxx b/gems/kvlGMMLikelihoodImageFilter.hxx
index b863226..daec95d 100644
--- a/gems/kvlGMMLikelihoodImageFilter.hxx
+++ b/gems/kvlGMMLikelihoodImageFilter.hxx
@@ -14,6 +14,11 @@ template< typename TInputImage >
 GMMLikelihoodImageFilter< TInputImage >
 ::GMMLikelihoodImageFilter()
 {
+#if ITK_VERSION_MAJOR >= 5
+  // use classic void ThreadedGenerateData( const OutputRegionType& threadRegion, ThreadIdType threadId )
+  // instead of default new signature void DynamicThreadedGenerateData( const OutputRegionType& threadRegion )
+  this->DynamicMultiThreadingOff();
+#endif  
 }
 
 
diff --git a/gems/kvlWMMLikelihoodImageFilter.hxx b/gems/kvlWMMLikelihoodImageFilter.hxx
index 8d24a44..c28b010 100644
--- a/gems/kvlWMMLikelihoodImageFilter.hxx
+++ b/gems/kvlWMMLikelihoodImageFilter.hxx
@@ -16,6 +16,11 @@ WMMLikelihoodImageFilter< TInputImage >
 ::WMMLikelihoodImageFilter()
 {
     //m_likelihoodFilter = GMMLikelihoodFilterType::New();
+#if ITK_VERSION_MAJOR >= 5 
+  // use classic void ThreadedGenerateData( const OutputRegionType& threadRegion, ThreadIdType threadId )
+  // instead of default new signature void DynamicThreadedGenerateData( const OutputRegionType& threadRegion )
+  this->DynamicMultiThreadingOff();
+#endif 
 }
 
 
diff --git a/samseg/cxx/module.cxx b/samseg/cxx/module.cxx
index 909d278..9e077cd 100644
--- a/samseg/cxx/module.cxx
+++ b/samseg/cxx/module.cxx
@@ -7,12 +7,15 @@
 #include "pyKvlTransform.h"
 #include "pyKvlRigidRegistration.h"
 #include "pyKvlAffineRegistration.h"
-#include "itkMultiThreader.h"
 
 namespace py = pybind11;
 
 void setGlobalDefaultNumberOfThreads(int maximumNumberOfThreads){
+#if ITK_VERSION_MAJOR >= 5
+    itk::MultiThreaderBase::SetGlobalDefaultNumberOfThreads( maximumNumberOfThreads );
+#else
     itk::MultiThreader::SetGlobalDefaultNumberOfThreads( maximumNumberOfThreads );
+#endif
 }
 
 PYBIND11_MODULE(gemsbindings, m) {
diff --git a/samseg/cxx/pyKvlCalculator.h b/samseg/cxx/pyKvlCalculator.h
index aaace71..645740d 100644
--- a/samseg/cxx/pyKvlCalculator.h
+++ b/samseg/cxx/pyKvlCalculator.h
@@ -155,7 +155,7 @@ class KvlCostAndGradientCalculator {
                 calculator->SetBoundaryCondition( kvl::AtlasMeshPositionCostAndGradientCalculator::SLIDING );
 
                 // Retrieve transform if one is provided
-                TransformType::ConstPointer  constTransform = 0;
+                TransformType::ConstPointer  constTransform = nullptr;
                 constTransform = static_cast< const TransformType* >( transform.m_transform.GetPointer() );
 
                 if ( constTransform.GetPointer() )