ransform parameters' gradient method of normalised cross-correlation

Algorithms
,
1

I found calculation transform parameters’ gradient method of normalised cross-correlation in itkNormalizedCorrelationImageToImageMetric.hxx
I don’t know why we calculate derivatives like this. I need you help to give me some suggestion

// Compute contributions to derivatives
  ti.GoToBegin();
  while (!ti.IsAtEnd())
  {
    index = ti.GetIndex();

    InputPointType inputPoint;
    fixedImage->TransformIndexToPhysicalPoint(index, inputPoint);

    if (this->m_FixedImageMask && !this->m_FixedImageMask->IsInsideInWorldSpace(inputPoint))
    {
      ++ti;
      continue;
    }

    OutputPointType transformedPoint = this->m_Transform->TransformPoint(inputPoint);

    if (this->m_MovingImageMask && !this->m_MovingImageMask->IsInsideInWorldSpace(transformedPoint))
    {
      ++ti;
      continue;
    }

    if (this->m_Interpolator->IsInsideBuffer(transformedPoint))
    {
      const RealType movingValue = this->m_Interpolator->Evaluate(transformedPoint);
      const RealType fixedValue = ti.Get();

      this->m_Transform->ComputeJacobianWithRespectToParametersCachedTemporaries(inputPoint, jacobian, jacobianCache);

      // Get the gradient by NearestNeighboorInterpolation:
      // which is equivalent to round up the point components.
      using CoordRepType = typename OutputPointType::CoordRepType;
      using MovingImageContinuousIndexType = ContinuousIndex<CoordRepType, MovingImageType::ImageDimension>;

      MovingImageContinuousIndexType tempIndex;
      this->m_MovingImage->TransformPhysicalPointToContinuousIndex(transformedPoint, tempIndex);

      typename MovingImageType::IndexType mappedIndex;
      mappedIndex.CopyWithRound(tempIndex);

      const GradientPixelType gradient = this->GetGradientImage()->GetPixel(mappedIndex);
      for (unsigned int par = 0; par < ParametersDimension; par++)
      {
        RealType sumF = NumericTraits<RealType>::ZeroValue();
        RealType sumM = NumericTraits<RealType>::ZeroValue();
        for (unsigned int dim = 0; dim < dimension; dim++)
        {
          const RealType differential = jacobian(dim, par) * gradient[dim];
          sumF += fixedValue * differential;
          sumM += movingValue * differential;
          if (this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0)
          {
            sumF -= differential * sf / this->m_NumberOfPixelsCounted;
            sumM -= differential * sm / this->m_NumberOfPixelsCounted;
          }
        }
        derivativeF[par] += sumF;
        derivativeM[par] += sumM;
      }
    }

    ++ti;
  }

  if (this->m_SubtractMean && this->m_NumberOfPixelsCounted > 0)
  {
    sff -= (sf * sf / this->m_NumberOfPixelsCounted);
    smm -= (sm * sm / this->m_NumberOfPixelsCounted);
    sfm -= (sf * sm / this->m_NumberOfPixelsCounted);
  }

  const RealType denom = -1.0 * std::sqrt(sff * smm);

  if (this->m_NumberOfPixelsCounted > 0 && denom != 0.0)
  {
    for (unsigned int i = 0; i < ParametersDimension; i++)
    {
      derivative[i] = (derivativeF[i] - (sfm / smm) * derivativeM[i]) / denom;
    }
  }
  else
  {
    for (unsigned int i = 0; i < ParametersDimension; i++)
    {
      derivative[i] = NumericTraits<MeasureType>::ZeroValue();
    }
  }
}