[Insight-users] Re: memory allocation for gradient

[Rupert Brooks]

Hi Lars,

I'm going to answer that in the general forum, in case others are
following along the conversation, and in case anyone spots errors in
what i say.

On 7/5/07, Lars Nygard <lnygard at yahoo.com> wrote:
[snipped]
> I now want to rewrite the MeansquareImageToImageMetric such that is uses
> the finite difference method. It seems to me that there is a difference where Image
> derivatives are being computed: the meansquaremetric has no place were image
> derivates are computed (this happen in the ImageToImageMetric in computeGradient),
> while the mattesmutualinformation does it using the finite difference method (in method
> ComputeImageDerivaties). Is this right? Where do you think I should insert the finite diff
> computation?

I did something like this - I apologize for not sharing my code
directly, but it isn't tested  perfectly yet, plus it is tied closely
to my research so theres a whole bunch of other stuff you would have
to take on... your better off to do it yourself.

Anyway, the following code excerpts might help a little, i hope.  I
approached it as a multi stage problem. So, i got each step working
before making the next change. In stage one, i just replaced the

in the mean squares metric with  compute image derivatives function
that does exactly the same thing.  But now, the gradient calculation
is isolated in one place, so i can change it at my convenience.
// replace ...

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

      typename MovingImageType::IndexType mappedIndex;
      for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ )
        {
        mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
        }

      const GradientPixelType gradient =
        this->GetGradientImage()->GetPixel( mappedIndex );

// with ...

	  GradientPixelType gradient;
      this->ComputeImageDerivatives(transformedPoint,gradient);


// and then elsewhere...

	template < class TFixedImage, class TMovingImage >
		void
		YourNewMetricImageToImageMetric2<TFixedImage,TMovingImage>
		::ComputeImageDerivatives(
		const MovingImagePointType& mappedPoint,
		GradientPixelType& gradient ) const
	{
		switch (m_GradientMethod) {
		case GradientImage:
			{
			// Get the gradient by NearestNeighboorInterpolation:
			// which is equivalent to round up the point components.
			typedef typename Superclass::OutputPointType OutputPointType;
			typedef typename OutputPointType::CoordRepType CoordRepType;
			typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
				MovingImageContinuousIndexType;

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

			typename MovingImageType::IndexType mappedIndex;
			for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ )
			{
				mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
			}
			gradient =
				this->GetGradientImage()->GetPixel( mappedIndex );
			}
			break;

This should give you numerically identical results to the old way.

Then - stage two - in the compute image derivatives function, i made
it switch on the type of derivative being computed and compute it the
Mattes info metric way in one case, old way in the other case.   Try
one way try the other, metric value shoud not change at all,
derivative will change small amounts due to numerical differences...
about 1% of the value or so.

// in the ComputeImageDerivatives method
	switch (m_GradientMethod) {
		case GradientImage:
			{
[.... as above...]
			break;

		case Calculator:
			{
			// For all generic interpolator use central differencing.
			gradient = m_DerivativeCalculator->Evaluate( mappedPoint );
			}
			break;
		

// dont forget to initialize the m_DerivativeCalculator in the
Initialize() method, and of course declare these variables and types
somewhere

// somewhere in the initialize method
	if(m_GradientMethod==Calculator) {
		m_DerivativeCalculator = DerivativeFunctionType::New();
		m_DerivativeCalculator->SetInputImage( m_MovingImage );
	}

Or heck, don't switch it, just change to using the calculator method.

So far that wont help your memory problem - you have to make the
Initialize method not compute the gradient image if the type of
gradient computation is not the regular kind.

// somewhere in the Initialize method.....
	if(m_GradientMethod==GradientImage) {
		m_ComputeGradient=true;
	} else {
		m_ComputeGradient=false;
	}

One last thing you might want to consider is to put these changes in a
base class and then inherit your NewMeanSquaresImageToImageMetric and
so on from it.  This is useful if you want to make this change for
other classes, such as the NormalizedCorrelationImageToImageMetric
which has the same problem.   Or just wait and eventually i will put
something in the Insight Journal.

Good luck,
Rupert B.