diff --git a/.github/workflows/build-test-package.yml b/.github/workflows/build-test-package.yml
index b4c70a334..ae2598b76 100644
--- a/.github/workflows/build-test-package.yml
+++ b/.github/workflows/build-test-package.yml
@@ -44,6 +44,7 @@ jobs:
         sudo apt-get install -y libgl1-mesa-glx
         sudo apt-get install -y libglvnd-core-dev
         sudo apt-get install -y freeglut3-dev
+        sudo apt-get install -y libtbb-dev
 
     - name: Install build dependencies
       run: |
@@ -201,7 +202,7 @@ jobs:
     strategy:
       max-parallel: 1
       matrix:
-        python-version: [36, 37, 38, 39]
+        python-version: [37, 38, 39]
         include:
           - itk-python-git-tag: "v5.3rc01"
             c-compiler: "gcc"
@@ -229,6 +230,7 @@ jobs:
         sudo apt-get install -y libgl1-mesa-glx
         sudo apt-get install -y libglvnd-core-dev
         sudo apt-get install -y freeglut3-dev
+        sudo apt-get install -y libtbb-dev
 
         python -m pip install --upgrade pip
         python -m pip install ninja
@@ -254,16 +256,16 @@ jobs:
       run: |
         cat > git-configure-build-sem.sh << EOF
         #!/bin/sh
-        ln -s -f -T /ITKPythonPackage/ITK-cp36-cp36m-manylinux2014_x64 /work/ITK-cp36-cp36m-manylinux2014_x64
+        ln -s -f -T /ITKPythonPackage/ITK-cp37-cp37m-manylinux2014_x86_64 /work/ITK-cp37-cp37m-manylinux2014_x86_64
         ln -s -f -T /ITKPythonPackage/ITK-source /work/ITK-source
-        export ITK_DIR="/work/ITK-cp36-cp36m-manylinux2014_x64"
+        export ITK_DIR="/work/ITK-cp37-cp37m-manylinux2014_x86_64"
         echo "itk_dir = ${ITK_DIR}"
         ls -la /work
         git clone https://github.com/Slicer/SlicerExecutionModel
         mkdir -p SlicerExecutionModel-build
         cd SlicerExecutionModel-build
         cmake -DCMAKE_BUILD_TYPE:STRING=Release \
-          -DITK_DIR=/work/ITK-cp36-cp36m-manylinux2014_x64 \
+          -DITK_DIR=/work/ITK-cp37-cp37m-manylinux2014_x86_64 \
           -GNinja \
           ../SlicerExecutionModel
         ninja
@@ -393,7 +395,7 @@ jobs:
     strategy:
       max-parallel: 1
       matrix:
-        python-version-minor: [6, 7, 8, 9]
+        python-version-minor: [7, 8, 9]
         include:
           - itk-python-git-tag: "v5.3rc01"
             c-compiler: "cl.exe"
diff --git a/src/Numerics/itktubeNJetImageFunction.hxx b/src/Numerics/itktubeNJetImageFunction.hxx
index f21372306..d9c4fb0c4 100644
--- a/src/Numerics/itktubeNJetImageFunction.hxx
+++ b/src/Numerics/itktubeNJetImageFunction.hxx
@@ -592,7 +592,7 @@ EvaluateAtContinuousIndex( const ContinuousIndexType & cIndex,
   // EVALUATE
   if( m_UseProjection )
     {
-    return EvaluateAtContinuousIndex( cIndex, scale );
+    m_MostRecentIntensity = EvaluateAtContinuousIndex( cIndex, scale );
     }
   else
     {
@@ -615,9 +615,8 @@ EvaluateAtContinuousIndex( const ContinuousIndexType & cIndex,
 
     m_MostRecentIntensity =  ( val0 + 0.5455*val1 + 0.5455*val2 )
       / ( 1+2*0.5455 ) / 2;
-
-    return m_MostRecentIntensity;
     }
+  return m_MostRecentIntensity;
 }
 
 
@@ -630,7 +629,7 @@ EvaluateAtContinuousIndex( const ContinuousIndexType & cIndex,
   // EVALUATE
   if( m_UseProjection )
     {
-    return EvaluateAtContinuousIndex( cIndex, scale );
+    m_MostRecentIntensity = EvaluateAtContinuousIndex( cIndex, scale );
     }
   else
     {
@@ -663,9 +662,8 @@ EvaluateAtContinuousIndex( const ContinuousIndexType & cIndex,
 
     m_MostRecentIntensity =
       ( 2*val0 + 0.5455*( val1+val2+val3+val4 ) )/( 2+4*0.5455 ) / 2;
-
-    return m_MostRecentIntensity;
     }
+  return m_MostRecentIntensity;
 }
 
 template< class TInputImage >
diff --git a/src/Registration/itktubeOptimizedSpatialObjectToImageRegistrationMethod.hxx b/src/Registration/itktubeOptimizedSpatialObjectToImageRegistrationMethod.hxx
index 9b6812847..badf420cb 100644
--- a/src/Registration/itktubeOptimizedSpatialObjectToImageRegistrationMethod.hxx
+++ b/src/Registration/itktubeOptimizedSpatialObjectToImageRegistrationMethod.hxx
@@ -175,7 +175,6 @@ void
 OptimizedSpatialObjectToImageRegistrationMethod<ObjectDimension, TImage>
 ::GenerateData( void )
 {
-  std::cout << "ORM: GenerateData" << std::endl;
   if( this->GetReportProgress() )
     {
     std::cout << "UPDATE START" << std::endl;
@@ -225,10 +224,8 @@ OptimizedSpatialObjectToImageRegistrationMethod<ObjectDimension, TImage>
       this->GetMovingSpatialObjectMaskObject() );
     }
 
-  std::cout << "ORM: GenerateData: Metric Init" << std::endl;
   metric->Initialize();
 
-  std::cout << "ORM: GenerateData: Optimize" << std::endl;
   try
     {
     this->Optimize(metric);
@@ -250,84 +247,98 @@ void
 OptimizedSpatialObjectToImageRegistrationMethod<ObjectDimension, TImage>
 ::Optimize( MetricType * metric )
 {
-  if( m_UseEvolutionaryOptimization && this->GetMaxIterations() > 0 )
+  if( this->GetMaxIterations() > 0 )
     {
-    if( this->GetReportProgress() )
-      {
-      std::cout << "EVOLUTIONARY START" << std::endl;
-      }
-
-    typedef OnePlusOneEvolutionaryOptimizer EvoOptimizerType;
-    EvoOptimizerType::Pointer evoOpt = EvoOptimizerType::New();
-
-    evoOpt->SetNormalVariateGenerator( Statistics::NormalVariateGenerator
-                                       ::New() );
-    evoOpt->SetEpsilon( this->GetTargetError() );
-    evoOpt->Initialize( 0.1 );
-    evoOpt->SetCatchGetValueException( true );
-    evoOpt->SetMetricWorstPossibleValue( 9999999 );
-    evoOpt->SetScales( this->GetTransformParametersScales() );
-    evoOpt->SetMaximumIteration( this->GetMaxIterations() );
-    evoOpt->SetMaximize( false );
-
-    if( this->GetObserver() )
-      {
-      evoOpt->AddObserver( IterationEvent(), this->GetObserver() );
-      }
-
-    if( this->GetReportProgress() )
+    if( m_UseEvolutionaryOptimization )
       {
-      typedef SpatialObjectToImageRegistrationViewer ViewerCommandType;
-      typename ViewerCommandType::Pointer command = ViewerCommandType::New();
-      if( this->GetTransform()->GetNumberOfParameters() > 16 )
+      if( this->GetReportProgress() )
         {
-        command->SetDontShowParameters( true );
+        std::cout << "EVOLUTIONARY START" << std::endl;
+        }
+  
+      typedef OnePlusOneEvolutionaryOptimizer EvoOptimizerType;
+      EvoOptimizerType::Pointer evoOpt = EvoOptimizerType::New();
+  
+      evoOpt->SetNormalVariateGenerator( Statistics::NormalVariateGenerator
+                                         ::New() );
+      evoOpt->SetEpsilon( this->GetTargetError() );
+      evoOpt->Initialize( 0.1 );
+      evoOpt->SetCatchGetValueException( true );
+      evoOpt->SetMetricWorstPossibleValue( 0 );
+      evoOpt->SetScales( this->GetTransformParametersScales() );
+      evoOpt->SetMaximumIteration( this->GetMaxIterations() * 0.75 );
+  
+      if( this->GetObserver() )
+        {
+        evoOpt->AddObserver( IterationEvent(), this->GetObserver() );
+        }
+  
+      if( this->GetReportProgress() )
+        {
+        typedef SpatialObjectToImageRegistrationViewer ViewerCommandType;
+        typename ViewerCommandType::Pointer command = ViewerCommandType::New();
+        if( this->GetTransform()->GetNumberOfParameters() > 16 )
+          {
+          command->SetDontShowParameters( true );
+          }
+        evoOpt->AddObserver( IterationEvent(), command );
+        }
+  
+      evoOpt->SetCostFunction(metric);
+      evoOpt->SetInitialPosition(m_InitialTransformParameters);
+      try
+        {
+        evoOpt->StartOptimization();
+        }
+      catch (...)
+        {
+        std::cout << "Exception caught during evolutionary registration."
+                  << std::endl;
+        std::cout << "Continuing using best values..." << std::endl;
+        std::cout << "  Pos = " << evoOpt->GetCurrentPosition()
+                  << std::endl << std::endl;
+        }
+  
+      if( evoOpt->GetCurrentPosition().size() ==
+        this->GetTransform()->GetNumberOfParameters() )
+        {
+        bool valid = true;
+        for(unsigned int i=0; i<evoOpt->GetCurrentPosition().size(); ++i)
+          {
+          if(evoOpt->GetCurrentPosition()[i] != evoOpt->GetCurrentPosition()[i])
+            {
+            valid = false;
+            break;
+            }
+          }
+        if( valid )
+          {
+          this->SetLastTransformParameters( evoOpt->GetCurrentPosition() );
+          this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
+          m_FinalMetricValue = metric->GetValue( m_LastTransformParameters );
+          }
+        else
+          {
+          std::cerr << "Error: Invalid Evolutionalry final parameters"
+            << std::endl;
+          this->SetLastTransformParameters( m_InitialTransformParameters );
+          this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
+          m_FinalMetricValue = 0;
+          }
+        }
+      if( this->GetReportProgress() )
+        {
+        std::cout << "EVOLUTIONARY END" << std::endl;
         }
-      evoOpt->AddObserver( IterationEvent(), command );
-      }
-
-    evoOpt->SetCostFunction(metric);
-    evoOpt->SetInitialPosition(m_InitialTransformParameters);
-    try
-      {
-      evoOpt->StartOptimization();
-      }
-    catch (...)
-      {
-      std::cout << "Exception caught during evolutionary registration."
-                << std::endl;
-      std::cout << "Continuing using best values..." << std::endl;
-      std::cout << "  Pos = " << evoOpt->GetCurrentPosition()
-                << std::endl << std::endl;
-      }
-
-    if( evoOpt->GetCurrentPosition().size() !=
-      this->GetTransform()->GetNumberOfParameters() )
-      {
-      this->SetLastTransformParameters( this->GetInitialTransformParameters() );
       }
     else
       {
-      this->SetLastTransformParameters( evoOpt->GetCurrentPosition() );
+      this->SetLastTransformParameters( m_InitialTransformParameters );
+      this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
+      m_FinalMetricValue = 0;
       }
-    this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
-    m_FinalMetricValue = metric->GetValue( m_LastTransformParameters );
-    }
-  else
-    {
-    this->SetLastTransformParameters( this->GetInitialTransformParameters() );
-    this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
-    m_FinalMetricValue = 0;
-    }
-
-  if( this->GetReportProgress() )
-    {
-    std::cout << "EVOLUTIONARY END" << std::endl;
-    }
 
-  if( this->GetMaxIterations() > 0 )
-    {
-    //if( this->GetReportProgress() )
+    if( this->GetReportProgress() )
       {
       std::cout << "GRADIENT START" << std::endl;
       }
@@ -335,15 +346,14 @@ OptimizedSpatialObjectToImageRegistrationMethod<ObjectDimension, TImage>
     typedef FRPROptimizer GradOptimizerType;
     GradOptimizerType::Pointer gradOpt = GradOptimizerType::New();
   
-    gradOpt->SetMaximize( false );
     gradOpt->SetCatchGetValueException( true );
-    gradOpt->SetMetricWorstPossibleValue( 9999999 );
-    gradOpt->SetStepLength( 0.25 );
+    gradOpt->SetMetricWorstPossibleValue( 0 );
+    gradOpt->SetStepLength( 0.1 );
     gradOpt->SetStepTolerance( this->GetTargetError() );
-    gradOpt->SetMaximumIteration( this->GetMaxIterations() );
+    gradOpt->SetMaximumIteration( this->GetMaxIterations() * 0.25 / 10 );
     gradOpt->SetMaximumLineIteration( 10 );
     gradOpt->SetScales( this->GetTransformParametersScales() );
-    gradOpt->SetUseUnitLengthGradient(true);
+    gradOpt->SetUseUnitLengthGradient(false);
     gradOpt->SetToFletchReeves();
   
     //if( this->GetReportProgress() )
@@ -362,7 +372,7 @@ OptimizedSpatialObjectToImageRegistrationMethod<ObjectDimension, TImage>
       }
   
     gradOpt->SetCostFunction(metric);
-    gradOpt->SetInitialPosition(m_InitialTransformParameters);
+    gradOpt->SetInitialPosition(m_LastTransformParameters);
     try
       {
       gradOpt->StartOptimization();
@@ -376,30 +386,42 @@ OptimizedSpatialObjectToImageRegistrationMethod<ObjectDimension, TImage>
                 << std::endl << std::endl;
       }
 
-    if( gradOpt->GetCurrentPosition().size() !=
+    if( gradOpt->GetCurrentPosition().size() ==
       this->GetTransform()->GetNumberOfParameters() )
       {
-      this->SetLastTransformParameters( this->GetInitialTransformParameters() );
+      bool valid = true;
+      for(unsigned int i=0; i<gradOpt->GetCurrentPosition().size(); ++i)
+        {
+        if(gradOpt->GetCurrentPosition()[i] != gradOpt->GetCurrentPosition()[i])
+          {
+          valid = false;
+          break;
+          }
+        }
+      if( valid )
+        {
+        this->SetLastTransformParameters( gradOpt->GetCurrentPosition() );
+        this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
+        m_FinalMetricValue = metric->GetValue( m_LastTransformParameters );
+        }
+      else
+        {
+        this->SetLastTransformParameters( m_InitialTransformParameters );
+        this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
+        m_FinalMetricValue = 0;
+        }
       }
-    else
+    if( this->GetReportProgress() )
       {
-      this->SetLastTransformParameters( gradOpt->GetCurrentPosition() );
+      std::cout << "GRADIENT END" << std::endl;
       }
-    this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
-    m_FinalMetricValue = metric->GetValue( m_LastTransformParameters );
     }
   else
     {
-    this->SetLastTransformParameters( this->GetInitialTransformParameters() );
+    this->SetLastTransformParameters( m_InitialTransformParameters );
     this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
     m_FinalMetricValue = 0;
     }
-
-
-  if( this->GetReportProgress() )
-    {
-    std::cout << "GRADIENT END" << std::endl;
-    }
 }
 
 template <unsigned int ObjectDimension, class TImage>
diff --git a/src/Registration/itktubePointBasedSpatialObjectToImageMetric.hxx b/src/Registration/itktubePointBasedSpatialObjectToImageMetric.hxx
index 896d4610f..89acfa2eb 100644
--- a/src/Registration/itktubePointBasedSpatialObjectToImageMetric.hxx
+++ b/src/Registration/itktubePointBasedSpatialObjectToImageMetric.hxx
@@ -561,9 +561,9 @@ PointBasedSpatialObjectToImageMetric< ObjectDimension, TFixedImage >
     value /= weightSum;
   }
 
-  std::cout << "GetValue : " << parameters << " = " << value << std::endl;
+  //std::cout << "GetValue : " << parameters << " = " << -value << std::endl;
 
-  return -value;
+  return -value;  // - to invert image so registration is minimization
 }
 
 template< unsigned int ObjectDimension, class TFixedImage >
@@ -653,6 +653,7 @@ PointBasedSpatialObjectToImageMetric< ObjectDimension, TFixedImage >
         tM = tM + outer_product( n2T, n2T );
         Vector<double, ImageDimension> tmpDeriv;
         imFunc->Derivative( fixedPoint, n1v, n2v, fixedScale, tmpDeriv );
+        value += *pointWeightIter * imFunc->Evaluate( fixedPoint, n1v, n2v, fixedScale );
         fixedDeriv.fill(0);
         for( unsigned int ii = 0; ii < ImageDimension; ++ii )
         {
@@ -664,6 +665,7 @@ PointBasedSpatialObjectToImageMetric< ObjectDimension, TFixedImage >
       {
         Vector<double, ImageDimension> tmpDeriv;
         imFunc->Derivative( fixedPoint, n1v, fixedScale, tmpDeriv );
+        value += *pointWeightIter * imFunc->Evaluate( fixedPoint, n1v, fixedScale );
         fixedDeriv.fill(0);
         for( unsigned int ii = 0; ii < ImageDimension; ++ii )
         {
@@ -675,7 +677,7 @@ PointBasedSpatialObjectToImageMetric< ObjectDimension, TFixedImage >
       tM = *pointWeightIter * tM;
       biasV += tM;
 
-      value += *pointWeightIter * imFunc->GetMostRecentIntensity();
+      //value += *pointWeightIter * imFunc->GetMostRecentIntensity();
 
       weightSum += *pointWeightIter;
     }
@@ -685,7 +687,7 @@ PointBasedSpatialObjectToImageMetric< ObjectDimension, TFixedImage >
 
   if( weightSum > 0 )
     {
-    value = - value / weightSum;
+    value = -value / weightSum;  // - to invert image
 
     VnlMatrixType biasVI = vnl_matrix_inverse< double >( biasV.as_ref() ).inverse();
 
@@ -709,7 +711,9 @@ PointBasedSpatialObjectToImageMetric< ObjectDimension, TFixedImage >
       {
         for( unsigned int d=0; d<ImageDimension; ++d)
         {
-          derivative[p] += jacobian[d][p] * -dXT[d];
+          derivative[p] -= jacobian[d][p] * dXT[d]; // - to invert image
+                                                    // so registration is
+                                                    // minimization.
         }
       }
   
@@ -717,6 +721,9 @@ PointBasedSpatialObjectToImageMetric< ObjectDimension, TFixedImage >
       ++fixedPointsDerivsIter;
     }
   }
+  //std::cout << "GetValueAndDerive : " << parameters << " = "
+  //          << -value << std::endl;
+  //std::cout << "                    " << derivative << std::endl;
 }
   
 template< unsigned int ObjectDimension, class TFixedImage >
@@ -783,6 +790,11 @@ bool
 PointBasedSpatialObjectToImageMetric< ObjectDimension, TFixedImage >
 ::IsValidFixedPoint( const FixedPointType & pnt ) const
 {
+  typename FixedImageType::IndexType indx;
+  if( !this->m_FixedImage->TransformPhysicalPointToIndex(pnt, indx) )
+    {
+    return false;
+    }
   if( this->m_FixedImageMaskObject.IsNotNull() &&
     this->m_UseFixedImageMaskObject )
     {
diff --git a/src/Registration/itktubeSpatialObjectToImageRegistrationHelper.hxx b/src/Registration/itktubeSpatialObjectToImageRegistrationHelper.hxx
index f19f4f798..a4eb561f7 100644
--- a/src/Registration/itktubeSpatialObjectToImageRegistrationHelper.hxx
+++ b/src/Registration/itktubeSpatialObjectToImageRegistrationHelper.hxx
@@ -478,7 +478,6 @@ SpatialObjectToImageRegistrationHelper<ObjectDimension, TImage>
     this->Initialize();
     }
 
-  std::cout << "HERE" << std::endl;
   if( m_EnableLoadedRegistration
       && m_LoadedMatrixTransform.IsNotNull() )
     {
@@ -499,7 +498,6 @@ SpatialObjectToImageRegistrationHelper<ObjectDimension, TImage>
     m_CurrentMatrixTransform = 0;
     }
 
-  std::cout << "HERE1" << std::endl;
   if( this->GetReportProgress() )
     {
     std::cout << "*** INITIAL REGISTRATION ***" << std::endl;
@@ -528,7 +526,6 @@ SpatialObjectToImageRegistrationHelper<ObjectDimension, TImage>
       regInit->SetMovingSpatialObjectMaskObject( m_MovingSpatialObjectMaskObject );
       }
     }
-  std::cout << "HERE2" << std::endl;
   if( m_EnableInitialRegistration )
     {
     switch( m_InitialMethodEnum )
@@ -551,7 +548,6 @@ SpatialObjectToImageRegistrationHelper<ObjectDimension, TImage>
       default:
         break;
       }
-    std::cout << "RH: Update: InitialReg" << std::endl;
     regInit->Update();
     m_InitialTransform = regInit->GetAffineTransform();
     }
@@ -570,7 +566,6 @@ SpatialObjectToImageRegistrationHelper<ObjectDimension, TImage>
       }
     }
 
-  std::cout << "RH: Update: Init done" << std::endl;
   m_CurrentMatrixTransform = m_InitialTransform;
 
   m_CompletedStage = INIT_STAGE;
@@ -578,7 +573,6 @@ SpatialObjectToImageRegistrationHelper<ObjectDimension, TImage>
 
   if( m_EnableRigidRegistration )
     {
-    std::cout << "RH: Update: registration" << std::endl;
     if( this->GetReportProgress() )
       {
       std::cout << "*** RIGID REGISTRATION ***" << std::endl;
@@ -642,7 +636,6 @@ SpatialObjectToImageRegistrationHelper<ObjectDimension, TImage>
       }
     regRigid->SetTransformParametersScales( scales );
 
-    std::cout << "RH: Update: set transform parameters scales" << std::endl;
     if( m_CurrentMatrixTransform.IsNotNull() )
       {
       regRigid->GetTypedTransform()->SetCenter(
@@ -657,9 +650,7 @@ SpatialObjectToImageRegistrationHelper<ObjectDimension, TImage>
         regRigid->GetTypedTransform()->GetParameters() );
       }
 
-    std::cout << "RH: Update: update" << std::endl;
     regRigid->Update();
-    std::cout << "RH: Update: update done" << std::endl;
 
     m_RigidTransform = RigidTransformType::New();
     m_RigidTransform->SetFixedParameters(
@@ -678,11 +669,8 @@ SpatialObjectToImageRegistrationHelper<ObjectDimension, TImage>
     m_CompletedResampling = false;
     }
 
-  std::cout << "RH: Update: done" << std::endl;
-
   if( m_EnableAffineRegistration )
     {
-    std::cout << "RH: Update: affine" << std::endl;
     this->AffineRegND<ImageDimension>();
     }
 }
@@ -764,13 +752,11 @@ SpatialObjectToImageRegistrationHelper<ObjectDimension, TImage>
         aTrans->GetParameters();
       typename MatrixTransformType::ParametersType tmpParams =
         tmpTrans->GetParameters();
-      std::cout << "Original params = " << aTransParams << std::endl;
       for( unsigned int p=0; p<tmpParams.size(); ++p )
         {
         tmpParams[p] = tmpParams[p] + portion * (aTransParams[p]-tmpParams[p]);
         }
       tmpTrans->SetParameters( tmpParams );
-      std::cout << "portion params = " << tmpParams << std::endl;
       }
     else
       {
@@ -874,7 +860,6 @@ SpatialObjectToImageRegistrationHelper<ObjectDimension, TImage>
   m_LoadedMatrixTransform->SetOffset( tfm.GetOffset() );
   if( invert )
     {
-    std::cout << "GetInverseTransform" << std::endl;
     typename MatrixTransformType::Pointer invTfm = MatrixTransformType::New();
     m_LoadedMatrixTransform->GetInverse( invTfm );
     m_LoadedMatrixTransform = invTfm;
diff --git a/test/Registration/CMakeLists.txt b/test/Registration/CMakeLists.txt
index 51e328974..bef0c7d86 100644
--- a/test/Registration/CMakeLists.txt
+++ b/test/Registration/CMakeLists.txt
@@ -165,7 +165,7 @@ itk_add_test(
     itktubeSpatialObjectToImageMetricTest
       DATA{${TubeTK_DATA_ROOT}/SyntheticVesselTubeImage.mha}
       DATA{${TubeTK_DATA_ROOT}/SyntheticVesselTubeManuallyModified.tre}
-      -0.0801461 )
+      -0.400731 )
 
 itk_add_test(
   NAME itktubeSpatialObjectToImageMetricTest2
diff --git a/test/Registration/itktubeSpatialObjectToImageRegistrationTest.cxx b/test/Registration/itktubeSpatialObjectToImageRegistrationTest.cxx
index 4a2a730e4..679c8ea94 100644
--- a/test/Registration/itktubeSpatialObjectToImageRegistrationTest.cxx
+++ b/test/Registration/itktubeSpatialObjectToImageRegistrationTest.cxx
@@ -67,7 +67,6 @@ int itktubeSpatialObjectToImageRegistrationTest( int argc, char * argv[] )
   typedef itk::tube::PointBasedSpatialObjectTransformFilter< TransformType, ObjectDimension >
                                                          TubeTransformFilterType;
 
-  std::cout << "start" << std::endl;
   // read image
   ImageReaderType::Pointer imageReader = ImageReaderType::New();
   imageReader->SetFileName( inputImage );
@@ -111,7 +110,6 @@ int itktubeSpatialObjectToImageRegistrationTest( int argc, char * argv[] )
     return EXIT_FAILURE;
     }
 
-  std::cout << "subsample" << std::endl;
   // subsample points in vessel
   typedef itk::tube::SubSampleSpatialObjectFilter< ObjectDimension > SubSampleFilterType;
   SubSampleFilterType::Pointer subSampleFilter =
@@ -174,12 +172,10 @@ int itktubeSpatialObjectToImageRegistrationTest( int argc, char * argv[] )
   std::cout << indent << translation[0] << " "
     << translation[1] << " " << translation[2] << std::endl;
 
-  double knownResult[] = { 0.00697,
-    -0.0049,
-    0.0052,
-    -1.8605,
-    -0.2859,
-    -0.9287 };
+  double knownResult[] = { 1.0, 0.0, 0.0,
+    0.0, 1.0, 0.0,
+    0.0, 0.0, 1.0
+    };
   int result = EXIT_SUCCESS;
   std::cout << "Parameters: " << std::endl;
   for( unsigned int ii = 0; ii < 6; ++ii )