[Insight-users] Memory concern
Luis Ibanez
luis . ibanez at kitware . com
Wed, 03 Sep 2003 17:42:38 -0400
Hi Stephan,
The large memory consumption of the metric is due to the internal use
of the GradientRecursiveGaussianImageFilter as an auxiliary for computing
the ImageGradient and with it the Metric derivatives.
This filters has N internal filters that allocate 'float' images.
(N=image dimension)
So the overhead per pixel of the input image is 4bytes * N. In
addition to that,
the resulting gradient image takes an image of N-D vectors in floats,
which is
4bytes * N. So by the time the gradient image is computed. The filter
has taken
2 * N * 4 bytes per pixel of input image.
In your case, an image of size 200x 200x200 pixels = 8x10^6 pixels
will result in an extra memory of : 2 * 3 * 4 * ( 8 x 10^6) = 192 x 10^6.
We could reduce the memory comsumption by clearing or destroying the
gradient filter, and holding a SmartPointer to the output vector image.
This will still leave 86x10^6 bytes of overhead.
You could also trade memory for performance and modify the MeanSquareMetric
for computing derivatives by using finite differences instead of using the
transformation Jacobian and the image gradient. This will probably result in
something like a factor 2X or 3X in speed.
Storing the gradient image as shorts instead of floats is also a possibility
and will let you win only 43x10^6 bytes for the gradient image. This still
will require to take (copy) the MeanSquares image metric and modify
the types used for storing the gradient image.
Note that the derivative of the image is actually computed by the metric
base class, the "ImageToImageMetric".
Yet another option is to (modify the metric) and keep a low resolution
version of the gradient image. The drawback of this option is that you
will be degrading the capabilities of the metric for guiding the
registration
towards a mimimum (optimum value).
----
As an experiment, I would suggest you to try the following:
1) disable the computation of the Gradient image in the
ImageToImageMetric.txx
(Insight/Code/Algorithms). In line 167 comment out the Update() call.
2) modify the MeanSquareMetrics (e.g. copy it with another name, both
for the file and the class inside). and just remove the GetDerivative
and GetValue an Derivative methods. When these two methods are not
overloaded, the ones in the base class are used. These base class
methods
compute derivatives by finite differences insteade of Jacobian *
ImageGradient.
3) Recompile your application and profile it again, both in speed and memory
consumption. You would probably find a factor 1/2X in memory and 3X to
4X in computing time.
Please let us know what you find. It that turns out to be reasonable
we may
want to add a public method to the ImageToImageMetric for disabling the
image gradient computations, and we should also consider adding a variant
of the MeanSquaresMetric to the set of ITK metrics. This variant will be
slower but less agressive in memory consumption.
In any case, if those changes get accepted, they will not be incorporated on
the toolkit until the release 1.4 is passed. The release is scheduled
for sept 15,
so for this eventual changes we should think in something like October.
When you perform the profiling of you application, watch carefully since
it is likely that the gradient filter is not the only responsible for those
317M....
Regards,
Luis
-------------------------
itk at stmoser . ch wrote:
>Hi
>
>I'm working within the same project as Michael. I have an application
>which does some evaluation of 3D registration runs that implements a
>function to calculate a metric value. If I work on my target data which
>is of type short and roughly 60M each image, my application allocates
>more than 800M of memory. I have set up a demo code that generates two
>200*200*200 images of type short and calculates a
>MeanSquaresImageToImageMetric with linear interpolation. Watching the
>application with 'top' I realize that it is 42M in size after creation
>of the images and rises up to 317M during the metric calculation.
>All the memory seems to be freed afterwards. Is there a way to cut peak
>memory usage?
>The actual problem is that depending on the hardware, the OS gets into
>swapping and makes calculations almost impossible.
>
>Thanks, Stephan
>
>------------------------------------------------------------------------
>
>#include "itkImage.h"
>#include "itkImageToImageMetric.h"
>#include "itkLinearInterpolateImageFunction.h"
>#include "itkMeanSquaresImageToImageMetric.h"
>#include "itkQuaternionRigidTransform.h"
>#include "itkBSplineInterpolateImageFunction.h"
>#include "itkVTKImageToImageFilter.h"
>#include "itkResampleImageFilter.h"
>#include "itkAffineTransform.h"
>
>typedef itk::Image<short, 3> Short3DImageType;
>typedef Short3DImageType::Pointer Short3DImageTypePointer;
>
>#define SIZE 200
>#define SPACING 1.0
>
>double CalculateMetricValue(Short3DImageTypePointer fixedImage,
> Short3DImageTypePointer movingImage)
>{
> typedef itk::ImageToImageMetric<Short3DImageType, Short3DImageType>
> ImageToImageMetricType;
> typedef ImageToImageMetricType::Pointer ImageToImageMetricTypePointer;
> typedef itk::QuaternionRigidTransform<double> TransformType;
>
> ImageToImageMetricTypePointer metric;
>
> typedef itk::MeanSquaresImageToImageMetric<Short3DImageType,
> Short3DImageType> MeanSquaresMetricType;
> metric = MeanSquaresMetricType::New();
>
> metric->SetFixedImage(fixedImage);
> metric->SetMovingImage(movingImage);
>
> TransformType::Pointer transform = TransformType::New();
> metric->SetTransform( transform );
>
> typedef itk::LinearInterpolateImageFunction<Short3DImageType,
> double> InterpolationType;
> InterpolationType::Pointer interpolator = InterpolationType::New();
> interpolator->SetInputImage( fixedImage );
> metric->SetInterpolator( interpolator );
>
> Short3DImageType::RegionType region =
> fixedImage->GetLargestPossibleRegion();
>
> metric->SetFixedImageRegion( region );
>
> ImageToImageMetricType::ParametersType params(7);
>
> params[0] = 0;
> params[1] = 0;
> params[2] = 1;
> params[3] = 1;
> params[4] = 0;
> params[5] = 0;
> params[6] = 0;
>
> metric->SetTransformParameters( params );
> try {
> metric->Initialize();
> metric->Modified();
> double metricValue = metric->GetValue( params );
>
> return metricValue;
> } catch (itk::ExceptionObject e) {
> cerr << e << endl;
> return 0;
> }
>}
>
>void GenerateData(Short3DImageTypePointer myImage) {
> const unsigned int Dimension = myImage->GetImageDimension();
> Short3DImageTypePointer::ObjectType::IndexType start;
> Short3DImageTypePointer::ObjectType::SizeType size;
> double* spacing = new double[Dimension];
> double* origin = new double[Dimension];
> for (int i = 0; i < Dimension; i++) {
> start[i] = 0;
> size[i] = SIZE;
> spacing[i] = SPACING;
> origin[i] = 0.0;
> }
> Short3DImageTypePointer::ObjectType::RegionType region;
> region.SetSize(size);
> region.SetIndex(start);
> myImage->SetRegions(region);
> myImage->Allocate();
> myImage->SetSpacing(spacing);
> myImage->SetOrigin(origin);
> Short3DImageTypePointer::ObjectType::IndexType idx;
> for(int x = 0; x < size[0]; x++) {
> for(int y = 0; y < size[1]; y++) {
> if (Dimension > 2) {
> // generate 3D data
> for(int z = 0; z < size[2]; z++) {
> idx[0] = x;
> idx[1] = y;
> idx[2] = z;
> if (x > size[0] / 4 && x < (size[0] * 3) / 4 &&
> y > size[1] / 4 && y < (size[1] * 3) / 4 &&
> z > size[2] / 4 && z < (size[2] * 3) / 4) {
> myImage->SetPixel(idx, 100);
> } else {
> myImage->SetPixel(idx, 0);
> }
> }
> } else {
> // generate 2D data
> idx[0] = x;
> idx[1] = y;
> if (x > size[0] / 4 && x < (size[0] * 3) / 4 &&
> y > size[1] / 4 && y < (size[1] * 3) / 4) {
> myImage->SetPixel(idx, 100);
> } else {
> myImage->SetPixel(idx, 0);
> }
> }
> }
> }
>}
>
>
>int main(int argc, char** argv) {
>
> Short3DImageTypePointer image1 = Short3DImageType::New();
> Short3DImageTypePointer image2 = Short3DImageType::New();
>
> cout << "setting up image 1..." << endl;
> GenerateData(image1);
> cout << "setting up image 2..." << endl;
> GenerateData(image2);
> cout << "calculating meansquaresmetric..." << endl;
> double metricValue = CalculateMetricValue(image1, image2);
>
> cout << metricValue << endl;
>
>}
>
>
>