ENH: Reorient gradients.

gandreadis · Mar 6, 2015 · 3a30614 · 3a30614
1 parent 63d9cc3
commit 3a30614
Show file tree

Hide file tree

Showing 3 changed files with 124 additions and 14 deletions.
diff --git a/Examples/itkantsRegistrationHelper.hxx b/Examples/itkantsRegistrationHelper.hxx
@@ -968,13 +968,21 @@ RegistrationHelper<TComputeType, VImageDimension>
 
           msqMetric->SetIntensityDistanceSigma( stageMetricList[currentMetricNumber].m_IntensityDistanceSigma );
           msqMetric->SetEuclideanDistanceSigma( stageMetricList[currentMetricNumber].m_EuclideanDistanceSigma );
-          if( msqMetric->GetEuclideanDistanceSigma() <= 0.0 || msqMetric->GetIntensityDistanceSigma() <= 0.0 )
+          if( msqMetric->GetEuclideanDistanceSigma() <= 0.0 )
             {
-            msqMetric->EstimateDistanceSigmasAutomaticallyOn();
+            msqMetric->EstimateEuclideanDistanceSigmaAutomaticallyOn();
             }
           else
             {
-            msqMetric->EstimateDistanceSigmasAutomaticallyOff();
+            msqMetric->EstimateEuclideanDistanceSigmaAutomaticallyOff();
+            }
+          if( msqMetric->GetIntensityDistanceSigma() <= 0.0 )
+            {
+            msqMetric->EstimateIntensityDistanceSigmaAutomaticallyOn();
+            }
+          else
+            {
+            msqMetric->EstimateIntensityDistanceSigmaAutomaticallyOff();
             }
 
           intensityPointSetMetric = msqMetric;

diff --git a/Temporary/itkMeanSquaresPointSetToPointSetIntensityMetricv4.h b/Temporary/itkMeanSquaresPointSetToPointSetIntensityMetricv4.h
@@ -72,10 +72,11 @@ class MeanSquaresPointSetToPointSetIntensityMetricv4:
   itkStaticConstMacro( FixedPointDimension, DimensionType, Superclass::FixedDimension );
 
   /**  Type of the moving point set. */
-  typedef TMovingPointSet                              MovingPointSetType;
-  typedef typename TMovingPointSet::PointType          MovingPointType;
-  typedef typename TMovingPointSet::PixelType          MovingPixelType;
-  typedef typename TMovingPointSet::PointsContainer    MovingPointsContainer;
+  typedef TMovingPointSet                                  MovingPointSetType;
+  typedef typename TMovingPointSet::PointType              MovingPointType;
+  typedef typename TMovingPointSet::PixelType              MovingPixelType;
+  typedef typename TMovingPointSet::PointsContainer        MovingPointsContainer;
+  typedef typename TMovingPointSet::PointDataContainer     MovingPointDataContainer;
 
   itkStaticConstMacro( MovingPointDimension, DimensionType, Superclass::MovingDimension );
 
@@ -96,13 +97,23 @@ class MeanSquaresPointSetToPointSetIntensityMetricv4:
   typedef typename Superclass::PointIdentifier         PointIdentifier;
   typedef typename Superclass::PointsConstIterator     PointsConstIterator;
 
+  typedef CovariantVector<TInternalComputationValueType, PointDimension> CovariantVectorType;
+
+  /**
+   * Set/get estimate intensity distance sigma automatically based on reasonable
+   * heuristics.
+   */
+  itkSetMacro( EstimateIntensityDistanceSigmaAutomatically, bool );
+  itkGetConstMacro( EstimateIntensityDistanceSigmaAutomatically, bool );
+  itkBooleanMacro( EstimateIntensityDistanceSigmaAutomatically );
+
   /**
-   * Set/get estimate distance sigmas automatically based on reasonable
+   * Set/get estimate Euclidean distance sigma automatically based on reasonable
    * heuristics.
    */
-  itkSetMacro( EstimateDistanceSigmasAutomatically, bool );
-  itkGetConstMacro( EstimateDistanceSigmasAutomatically, bool );
-  itkBooleanMacro( EstimateDistanceSigmasAutomatically );
+  itkSetMacro( EstimateEuclideanDistanceSigmaAutomatically, bool );
+  itkGetConstMacro( EstimateEuclideanDistanceSigmaAutomatically, bool );
+  itkBooleanMacro( EstimateEuclideanDistanceSigmaAutomatically );
 
   /**
    * Set/get intensity sigma -- modulate the intensity distance contribution in the
@@ -146,17 +157,36 @@ class MeanSquaresPointSetToPointSetIntensityMetricv4:
   MeanSquaresPointSetToPointSetIntensityMetricv4();
   virtual ~MeanSquaresPointSetToPointSetIntensityMetricv4();
 
+  /**
+   * Estimate the intensity distance sigma based on simple heuristic
+   */
   void EstimateIntensityDistanceSigma();
+
+  /**
+   * Estimate the Euclidean distance sigma based on simple heuristic
+   */
   void EstimateEuclideanDistanceSigma();
 
+  /**
+   * Warp the fixed point set gradients based on the fixed transform.
+   */
+  void TransformFixedPointSetGradients() const;
+
+  /**
+   * Warp the moving point set gradients based on the moving transform.
+   */
+  void TransformMovingPointSetGradients() const;
+
+
   /** PrintSelf function */
   void PrintSelf( std::ostream & os, Indent indent ) const ITK_OVERRIDE;
 
 private:
   MeanSquaresPointSetToPointSetIntensityMetricv4(const Self &); //purposely not implemented
   void operator=(const Self &);               //purposely not implemented
 
-  bool                                m_EstimateDistanceSigmasAutomatically;
+  bool                                m_EstimateIntensityDistanceSigmaAutomatically;
+  bool                                m_EstimateEuclideanDistanceSigmaAutomatically;
   TInternalComputationValueType       m_IntensityDistanceSigma;
   TInternalComputationValueType       m_EuclideanDistanceSigma;
 

diff --git a/Temporary/itkMeanSquaresPointSetToPointSetIntensityMetricv4.hxx b/Temporary/itkMeanSquaresPointSetToPointSetIntensityMetricv4.hxx
@@ -31,7 +31,8 @@ MeanSquaresPointSetToPointSetIntensityMetricv4<TFixedPointSet, TMovingPointSet,
   this->m_EuclideanDistanceSigma = std::sqrt( 5.0 );
   this->m_IntensityDistanceSigma = std::sqrt( 5.0 );
 
-  this->m_EstimateDistanceSigmasAutomatically = true;
+  this->m_EstimateIntensityDistanceSigmaAutomatically = true;
+  this->m_EstimateEuclideanDistanceSigmaAutomatically = true;
 
   this->m_UsePointSetData = true;
 }
@@ -51,13 +52,84 @@ MeanSquaresPointSetToPointSetIntensityMetricv4<TFixedPointSet, TMovingPointSet,
 {
   Superclass::Initialize();
 
-  if( this->m_EstimateDistanceSigmasAutomatically )
+  this->TransformMovingPointSetGradients();
+  this->TransformFixedPointSetGradients();
+
+  if( this->m_EstimateIntensityDistanceSigmaAutomatically )
     {
     this->EstimateIntensityDistanceSigma();
+    }
+  if( this->m_EstimateEuclideanDistanceSigmaAutomatically )
+    {
     this->EstimateEuclideanDistanceSigma();
     }
 }
 
+template<typename TFixedPointSet, typename TMovingPointSet, class TInternalComputationValueType>
+void
+MeanSquaresPointSetToPointSetIntensityMetricv4<TFixedPointSet, TMovingPointSet, TInternalComputationValueType>
+::TransformFixedPointSetGradients() const
+{
+  // Transform the moving point set data with the moving transform.
+  // We calculate the value and derivatives in the moving space.
+
+  typename FixedPointsContainer::ConstIterator It = this->m_FixedPointSet->GetPoints()->Begin();
+
+  while( It != this->m_FixedPointSet->GetPoints()->End() )
+    {
+    PixelType pixel;
+    NumericTraits<PixelType>::SetLength( pixel, 1 );
+    bool doesPointDataExist = this->m_FixedPointSet->GetPointData( It.Index(), &pixel );
+    if( ! doesPointDataExist )
+      {
+      itkExceptionMacro( "The corresponding data for point " << It.Value() << " (pointId = " << It.Index() << ") does not exist." );
+      }
+    SizeValueType numberOfVoxelsInNeighborhood = pixel.size() / ( 1 + PointDimension );
+
+    for( SizeValueType n = 0; n < numberOfVoxelsInNeighborhood; n++ )
+      {
+      CovariantVectorType covariantVector;
+
+      for( unsigned int d = 0; d < PointDimension; d++ )
+        {
+        covariantVector[d] = pixel[n * ( PointDimension + 1 ) + d + 1];
+        }
+
+      // Here we assume that transforming the vector at the neighborhood voxel
+      // is close to performing the transformation at the center voxel.
+
+      CovariantVectorType transformedCovariantVector =
+        this->m_FixedTransform->TransformCovariantVector( covariantVector, It.Value() );
+
+      for( unsigned int d = 0; d < PointDimension; d++ )
+        {
+        pixel[n * ( PointDimension + 1 ) + d + 1] = transformedCovariantVector[d];
+        }
+      }
+
+    // evaluation is perfomed in moving space, so just copy
+    this->m_FixedTransformedPointSet->SetPointData( It.Index(), pixel );
+    ++It;
+    }
+}
+
+template<typename TFixedPointSet, typename TMovingPointSet, class TInternalComputationValueType>
+void
+MeanSquaresPointSetToPointSetIntensityMetricv4<TFixedPointSet, TMovingPointSet, TInternalComputationValueType>
+::TransformMovingPointSetGradients() const
+{
+  // Transform the moving point set data with the moving transform.
+  // We calculate the value and derivatives in the moving space.
+
+  typename MovingPointDataContainer::ConstIterator It = this->m_MovingPointSet->GetPointData()->Begin();
+  while( It != this->m_MovingPointSet->GetPointData()->End() )
+    {
+    // evaluation is perfomed in moving space, so just copy
+    this->m_MovingTransformedPointSet->SetPointData( It.Index(), It.Value() );
+    ++It;
+    }
+}
+
 template<typename TFixedPointSet, typename TMovingPointSet, class TInternalComputationValueType>
 void
 MeanSquaresPointSetToPointSetIntensityMetricv4<TFixedPointSet, TMovingPointSet, TInternalComputationValueType>