Neighborhood Cross Correlation

Algorithms
1

This is a function in itkANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader.hxx to compute cross correlation in a local window. I think if fabs(sFixedFixed_sMovingMoving) < NumericTraits< LocalRealType >::epsilon(), localCC should not be one (which means perfect correlation). For example, [0, 0, 0] and [1, 2, 3] are actually not correlated, but this function will give localCC the value of one.

template < typename TDomainPartitioner, typename TImageToImageMetric, typename TNeighborhoodCorrelationMetric >
void
ANTSNeighborhoodCorrelationImageToImageMetricv4GetValueAndDerivativeThreader< TDomainPartitioner, TImageToImageMetric, TNeighborhoodCorrelationMetric >
::ComputeMovingTransformDerivative( const ScanIteratorType &, ScanMemType &scanMem, const ScanParametersType &, DerivativeType &deriv, MeasureType &localCC, const ThreadIdType threadId) const
{
MovingImageGradientType derivWRTImage;
localCC = NumericTraits::OneValue();

typedef InternalComputationValueType LocalRealType;

LocalRealType sFixedFixed = scanMem.sFixedFixed;
LocalRealType sMovingMoving = scanMem.sMovingMoving;
LocalRealType sFixedMoving = scanMem.sFixedMoving;
LocalRealType fixedI = scanMem.fixedA;
LocalRealType movingI = scanMem.movingA;

LocalRealType sFixedFixed_sMovingMoving = sFixedFixed * sMovingMoving;

if ( fabs(sFixedFixed_sMovingMoving) > NumericTraits< LocalRealType >::epsilon() )
{
localCC = sFixedMoving * sFixedMoving / (sFixedFixed_sMovingMoving);
}

if( this->m_ANTSAssociate->GetComputeDerivative() )
{
const MovingImageGradientType movingImageGradient = scanMem.movingImageGradient;

if ( ! (sFixedFixed > NumericTraits<LocalRealType>::epsilon() && sMovingMoving > NumericTraits<LocalRealType>::epsilon() ) )
  {
  deriv.Fill( NumericTraits<DerivativeValueType>::ZeroValue() );
  return;
  }

for (ImageDimensionType qq = 0; qq < TImageToImageMetric::VirtualImageDimension; qq++)
  {
  derivWRTImage[qq] = 2.0 * sFixedMoving / (sFixedFixed_sMovingMoving) * (fixedI - sFixedMoving / sMovingMoving * movingI) * movingImageGradient[qq];
  }

/* Use a pre-allocated jacobian object for efficiency */
typedef JacobianType & JacobianReferenceType;
JacobianReferenceType jacobian = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobian;
JacobianReferenceType jacobianPositional = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobianPositional;

/** For dense transforms, this returns identity */
this->m_Associate->GetMovingTransform()->
  ComputeJacobianWithRespectToParametersCachedTemporaries(scanMem.virtualPoint,
                                                          jacobian,
                                                          jacobianPositional);

NumberOfParametersType numberOfLocalParameters = this->m_Associate->GetMovingTransform()->GetNumberOfLocalParameters();

for (NumberOfParametersType par = 0; par < numberOfLocalParameters; par++)
  {
  deriv[par] = NumericTraits<DerivativeValueType>::ZeroValue();
  for (ImageDimensionType dim = 0; dim < TImageToImageMetric::MovingImageDimension; dim++)
    {
    deriv[par] += derivWRTImage[dim] * jacobian(dim, par);
    }
  }
}

}