[Insight-users] registration and spacing: GetOutput() / SetSpacing()
and const-correctness
Luis Ibanez
luis . ibanez at kitware . com
Mon, 14 Jul 2003 10:27:21 -0400
Hi Michael,
Thanks for sending your code.
It helps to clarify the produre you are using.
Here are some factor that can be preventing your
registration process from working correctly.
1) The mechanism you are using for changing the
image scale is not reliable. You are taking
the output of the reader and invoking its
SetSpacing() method.
double spacing[2];
spacing[0] = 1.0 * FACTOR;
spacing[1] = 1.0 * FACTOR;
reader->GetOutput()->SetSpacing(spacing);
*img = reader->GetOutput();
You shouldn't attempt to change the spacing
of the Output() in any ITK filter, nor any
other property of the filter output. The output
of a filter belongs to the filter itself, and
can be recomputed at any call of the Update()
method. The changes of spacing that you
are applying could be overriden by the next
call to an Update() in the pipeline.
If you double check in your GetAIM() method,
you are invoking Update() on the importer,
after having reset the Spacing of the reader.
This Update() will override your spacing
modifications.
The GetOutput() method of ITK filters should
infact return a ConstSmartPointer in order
to prevent this miss-use of the output of a filter.
This is something we may have to fix in a future
release...
Since you are doing this spacing change more
for debugging purposes than for being part of
a normal data pipeline, I'll suggest you to
convert your PNG images into MetaImage format
and then simply modify the spacing in the .mhd
file. Then, execute your program using the .mhd
files as input (instead of the PNG).
Another option is to use the "ChangeInformationImageFilter"
which allows you to modify image data in a pipeline
compatible fashion. Note that this is a very risky
operation that you would not like to have in a
real medical application. Please don't use this
for anything different than debugging.
2) In the conversion from VTK image to ITK image
you may find convenient to use the filters:
ImageToVTKImageFilter
VTKImateToImageFilter
available in InsightApplications/Auxiliary/vtk
That will simplify the code in your GetAIM() function.
3) The only occasion in which you are entitled to call
the SetSpacing() and SetOrigin() methods is when you
create the image in a program, instead of reading it
from a file, or taking it from the output of a filter.
Filter outputs should not be modified.
4) PNG is a poor format for storing medical images.
It is supported in ITK only for facilitating simple
experiments, examples and demos. In real applications
you may want to use a more reliable format like GIPL,
Analyze, MetaImage or VTK.
Regards,
Luis
------------------------
Michael Kuhn wrote:
> Hi Luis,
>
> I think, my explanation of the setup was not very clear. I've attached
> my code to avoid any misunderstandings.
>
> I do in fact scale the spacing of both images, the fixed as well as the
> moving image. The two final Print statements in the code show the same
> spacing for any value of FACTOR.
>
> We introduced this FACTOR, since we noticed that the registration worked
> for images with spacing 1.0 but didn't for images with spacing 0.034.
> Therefore we tried to alter the working version (using spacing 1.0) into
> a version with different spacing and expected it to work. However, it
> didn't.
>
> Since transformations are done in physical space, we scale the transform
> vector as well as the initial parameters and the step sizes. Is there
> something more we would have to scale by the same factor?
>
> In case the program is fed with a 2D png image, I write the fixed and
> the moving images to disk. The resulting images seem to be ok for me
> (i.e. one looks exactly like the original one and the other has the same
> dimensions (in pixels) but is translated a little bit to the top right
> and is black at the left and bottom edge).
>
> In case of 2D png images, the registrations produces incorrect results
> for FACTOR above approximately 4 and it doesn't iterate (but generates
> an end event immediately after the start event) for FACTOR values below
> apprixunately 0.1.
>
> In case of 3D aim data (an image type used for CT images by some people,
> which uses short intensities), the behaviour is slightly different. In
> this case, the registrations always iterates, but produces wrong results
> (except for spacing 1.0).
>
> Thanks a lot for your help,
>
> Michael
>
>
>
> Luis Ibanez wrote:
>
>>
>> Hi Michael,
>>
>> It is normal that the registration fails in the setup that you
>> are configuring.
>>
>> Your ImageRegistrationMethod is using a "TranslationTransform".
>>
>> This transform cannot compensate for the scaling changes that
>> you are introducing when you change the pixel spacing.
>>
>> If you want to compensate for translation and scaling you may have
>> to use a "SimilarityTransform" which is the simplest transform
>> capable of managing both translation and scaling. You could also
>> go one step further and use an "AffineTransform".
>>
>> Your expectation of the registration to work with a common
>> scaling of the factors you altered, seems to be based on the
>> assumption that registraion is performed in pixel space. That
>> is, in the discrete grid of the images.
>>
>> However, registration in ITK is not performed in the pixel space.
>> It is performed in the physical space in millimeters.
>>
>> When you scale the pixel spacing of one image, you create the
>> conditions for registering two objects that have different
>> physical sizes. The correspondance between these two object
>> can only be achieved with a Transform supporting some form of
>> scaling factor.
>>
>>
>> ---
>>
>>
>> Just on the side: Note that a scale factor of 2X is quite
>> aggresive. In normal conditions in medical image registration
>> we would not expect to find such a huge scale factor in normal
>> anatomical structures, unless you were comparing infants with
>> adults. E.g. For reference: In humans, the size of the brain
>> in the newborn is about 25% of its adult size and it reaches
>> about 75 % of its adult size by the end of the first year.
>>
>>
>>
>>
>> Regards
>>
>>
>>
>> Luis
>>
>>
>> -----------------------
>> Michael Kuhn wrote:
>>
>>> Hi,
>>>
>>> I'm trying to do a registration using 3D Data (pixel type short). The
>>> general setup is as follows:
>>>
>>> The input consists of one file containing the data. The data is then
>>> transformed using a ResampleImageFilter which applies a translation
>>> to it. Afterwards, the original as well as the transformed image are
>>> fed into a ImageRegistrationMethod as fixed (original) and moving
>>> (translated) image respectively.
>>>
>>> If the spacing of both data sets is set to 1.0 (in every direction),
>>> the registration works fine. However if I scale the spacings as well
>>> as the translation vector (for the translation itself as well as the
>>> initalParameters for the registration) and the minimum/maximum step
>>> length with the same factor, the registration doesn't work any longer
>>> (i.e. produces incorrect results). The origin of the data was set to
>>> [0,0,0].
>>>
>>> I would expect that the registration works for any common scaling of
>>> these parameters, however, it doesn't. Can somebody explain me why?
>>>
>>> My dataset has a size of [40,40,40] and I translate it with
>>> parameters FACTOR * [-1, -2, 0]. The initial parameters are set to
>>> FACTOR * [2,2,0] (where FACTOR is the common scaling factor mentioned
>>> above). The max and min step length are set to FACTOR * 4.0 and
>>> FACTOR * 0.01 respectively.
>>>
>>> For the initial transform as well as the registration I use a
>>> TranslateTransform<double, 3> and a
>>> LinearInterpolateImageFunction<ImageType, double>. The registration
>>> is done using a RegularStepGradientDescentOptimzer and a
>>> MeanSquaresImageToImageMetric.
>>>
>>> Thanks,
>>>
>>> Michael
>>>
>>> _______________________________________________
>>> Insight-users mailing list
>>> Insight-users at itk . org
>>> http://www . itk . org/mailman/listinfo/insight-users
>>>
>>
>>
>>
>>
>
>
> ------------------------------------------------------------------------
>
> #include "itkRGBPixel.h"
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
> #include "itkRescaleIntensityImageFilter.h"
> #include "itkImage.h"
>
> #include "itkImageRegistrationMethod.h"
> #include "itkTranslationTransform.h"
> #include "itkMeanSquaresImageToImageMetric.h"
> #include "itkLinearInterpolateImageFunction.h"
> #include "itkBSplineInterpolateImageFunction.h"
> #include "itkRegularStepGradientDescentOptimizer.h"
> #include "itkNearestNeighborInterpolateImageFunction.h"
>
> #include "itkResampleImageFilter.h"
>
> #include "itkVTKImageImport.h"
>
> #include "vtkAIMReader.h"
> #include "vtkImageExport.h"
>
> #include "vtkUtility.h"
> #include "vtkPipelineConnectorVTK2ITK.h"
> #include "observer.h"
>
> #include <iostream>
>
> #define FACTOR 2.0
>
> typedef short PixelType;
> typedef itk::Image<PixelType, 2> ImageType;
>
> template<typename ImageTypePointer>
> void GetPNG(char* inputFilename, ImageTypePointer *img) {
> typedef itk::ImageFileReader<ImageTypePointer::ObjectType> ReaderType;
> ReaderType::Pointer reader = ReaderType::New();
> reader->SetFileName(inputFilename);
> reader->Update();
> double spacing[2];
> spacing[0] = 1.0 * FACTOR;
> spacing[1] = 1.0 * FACTOR;
> reader->GetOutput()->SetSpacing(spacing);
> *img = reader->GetOutput();
> }
>
> template<typename ImageTypePointer>
> void GetAIM(char* inputFilename, ImageTypePointer *img) {
> vtkAIMReader *reader = vtkAIMReader::New();
> reader->SetFileName(inputFilename);
> reader->Update();
> float* spacing = reader->GetOutput()->GetSpacing();
> cout << "aimreader.spacing: " << spacing[0] << " " << spacing[1] <<
> " " << spacing[2] << endl;
> spacing[0] = FACTOR * 1.0; //FACTOR * 0.019998;
> spacing[1] = FACTOR * 1.0; //FACTOR * 0.019998;
> spacing[2] = FACTOR * 1.0; //FACTOR * 0.021;
> reader->GetOutput()->SetSpacing(spacing);
> float* origin = reader->GetOutput()->GetOrigin();
> origin[0] = FACTOR * 0.0;
> origin[1] = FACTOR * 0.0;
> origin[2] = FACTOR * 0.0;
> reader->GetOutput()->SetOrigin(origin);
> reader->GetOutput()->PrintSelf(std::cout, 0);
>
> // reader->Update();
> vtkImageExport *vtkExporter = vtkImageExport::New();
> vtkExporter->SetInput(reader->GetOutput());
> typedef itk::VTKImageImport<ImageTypePointer::ObjectType>
ImageImportType;
> ImageImportType::Pointer itkImporter = ImageImportType::New();
> vtkPipelineConnectorVTK2ITK<vtkImageExport*, ImageImportType>::
> ConnectPipelines(vtkExporter, itkImporter);
> itkImporter->Update();
> *img = itkImporter->GetOutput();
> }
>
> template <typename ImageTypePointer>
> void GenerateData(ImageTypePointer* myImage) {
> const unsigned int Dimension = (*myImage)->GetImageDimension();
> ImageTypePointer::ObjectType::IndexType start;
> ImageTypePointer::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] = 40;
> spacing[i] = 1.0;
> origin[i] = 0.0;
> }
> ImageTypePointer::ObjectType::RegionType region;
> region.SetSize(size);
> region.SetIndex(start);
> (*myImage)->SetRegions(region);
> (*myImage)->Allocate();
> (*myImage)->SetSpacing(spacing);
> (*myImage)->SetOrigin(origin);
> ImageTypePointer::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);
> }
> }
> }
> }
> } // generate own data
>
> char* GetInterpolatorType(int argc, char** argv) {
> char* type = new char[20];
> if (vtkUtility::GetParam("-i", type, argc, argv)) {
> if (strcmp(type, "linear") == 0 || strcmp(type, "bspline") == 0 ||
> strcmp(type, "nearest") == 0) {
> return type;
> }
> }
> return "linear";
> }
>
> template <typename RegistrationTypePointer>
> void SetInterpolator(RegistrationTypePointer* reg, char* type) {
> if (strcmp(type, "linear") == 0) {
> typedef itk::LinearInterpolateImageFunction<
> RegistrationTypePointer::ObjectType::MovingImageType, double>
> InterpolatorType;
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> (*reg)->SetInterpolator(interpolator);
> cout << "using linear interpolator." << endl;
> } else if (strcmp(type, "bspline") == 0) {
> // http:// compare: public.kitware.com/Insight/Web/Doxygen/html/
> // classitk_1_1BSplineInterpolateImageFunction.html#s10
> // for the templating
> typedef itk::BSplineInterpolateImageFunction<
> RegistrationTypePointer::ObjectType::MovingImageType, double,
> double> InterpolatorType;
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> (*reg)->SetInterpolator(interpolator);
> interpolator->SetSplineOrder(3);
> cout << "using bspline interpolator with order 3." << endl;
> } else if (strcmp(type, "nearest") == 0) {
> typedef itk::NearestNeighborInterpolateImageFunction<
> RegistrationTypePointer::ObjectType::MovingImageType, double>
> InterpolatorType;
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> (*reg)->SetInterpolator(interpolator);
> cout << "using nearest neighbor interpolator." << endl;
> }
> }
>
>
> template <typename ImageTypePointer>
> itk::InterpolateImageFunction::Pointer
GetInterpolator(ImageTypePointer img) {
> if (strcmp(typ, "bspline") == 0) {
> // http:// compare: public.kitware.com/Insight/Web/Doxygen/html/
> // classitk_1_1BSplineInterpolateImageFunction.html#s10
> // for the templating
> typedef itk::BSplineInterpolateImageFunction<
> ImageTypePointer::ObjectType, double, double>
> InterpolatorType;
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> cout << "using bspline interpolator with order 3." << endl;
> interpolator->SetSplineOrder(3);
> return interpolator;
> } else if (strcmp(type, "nearest") == 0) {
> typedef itk::NearestNeighborInterpolateImageFunction<
> ImageTypePointer::ObjectType, double>
> InterpolatorType;
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> cout << "using nearest neighbor interpolator." << endl;
> return interpolator;
> }
> // default: linear interpolator
> typedef
itk::LinearInterpolateImageFunction<ImageTypePointer::ObjecType,
> double>
> InterpolatorType;
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> cout << "using linear interpolator." << endl;
> return interpolator;
> }
>
>
> int main(int argc, char** argv)
> {
> char inputFilename[255];
>
> const unsigned int Dimension = 2; // use 2 in case of png images
> // use 3 in case of aim images
> // typedef short PixelType;
> // typedef itk::Image<PixelType, Dimension> ImageType;
>
> bool png = false;
>
> ImageType::Pointer myImage;
> if (vtkUtility::GetParam("-png", inputFilename, argc, argv)) {
> GetPNG(inputFilename, &myImage);
> cout << "PNG image imported." << endl;
> png = true;
> } else if (vtkUtility::GetParam("-aim", inputFilename, argc, argv)) {
> GetAIM(inputFilename, &myImage);
> cout << "aim image imported" << endl;
> } else {
> //generate own data
> myImage = ImageType::New();
> GenerateData(&myImage);
> cout << "own data generated" << endl;
> }
>
> myImage->Print(std::cout);
>
> typedef itk::ResampleImageFilter<ImageType, ImageType> FilterType;
> FilterType::Pointer filter = FilterType::New();
>
>
>
> typedef itk::TranslationTransform<double, Dimension> TransformType;
> typedef itk::LinearInterpolateImageFunction<ImageType, double>
> InterpolatorType;
> InterpolatorType::Pointer myInterpolator = InterpolatorType::New();
> filter->SetInterpolator(myInterpolator);
>
>
> filter->SetDefaultPixelValue(0);
>
> //double spacing[Dimension];
> //double origin[Dimension];
> //for (int i = 0; i < Dimension; i++) {
> // spacing[i] = 1.0;
> // origin[i] = 1.0;
> //}
>
> const double* spacing = myImage->GetSpacing();
> const double* origin = myImage->GetOrigin();
>
> filter->SetOutputSpacing(spacing);
> filter->SetOutputOrigin(origin);
>
> ImageType::SizeType size =
> myImage->GetLargestPossibleRegion().GetSize();
> //size[0] = size[0] + 10;
> //size[1] = size[1] + 15;
>
> std::cout << "size: " << size[0] << " " << size[1] << std::endl;
> filter->SetSize(size);
>
> TransformType::Pointer myTransform = TransformType::New();
> TransformType::OutputVectorType translation;
> translation[0] = FACTOR * -1;
> translation[1] = FACTOR * -2;
> if (Dimension > 2) {
> translation[2] = 0.0;
> }
>
> myTransform->Translate(translation);
> filter->SetTransform(myTransform);
>
> filter->SetInput(myImage);
>
> typedef itk::ImageRegistrationMethod<ImageType, ImageType>
> RegistrationType;
> RegistrationType::Pointer registration = RegistrationType::New();
>
> //****************************************
>
> typedef itk::LinearInterpolateImageFunction<ImageType, double>
> InterpolatorType;
>
> // InterpolatorType::Pointer interpolator =
GetInterpolator(registration->GetMovingImage());
>
> //****************************************
>
>
> SetInterpolator(®istration, GetInterpolatorType(argc, argv));
> char* type = GetInterpolatorType(argc, argv);
> // registration->SetInterpolator(interpolator);
>
> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
> typedef itk::MeanSquaresImageToImageMetric<ImageType, ImageType>
> MetricType;
>
> MetricType::Pointer metric = MetricType::New();
> TransformType::Pointer transform = TransformType::New();
> OptimizerType::Pointer optimizer = OptimizerType::New();
> //InterpolatorType::Pointer interpolator = InterpolatorType::New();
>
> registration->SetMetric(metric);
> registration->SetOptimizer(optimizer);
> registration->SetTransform(transform);
> //registration->SetInterpolator(interpolator);
>
> registration->SetFixedImage(myImage);
>
> cout << "%%%%%%%%%%%%%%fixed image%%%%%%%%%%%%%%%%%" << endl;
> myImage->Print(std::cout);
> cout << endl;
>
> filter->Update();
>
> registration->SetMovingImage(filter->GetOutput());
> cout << "%%%%%%%%%%%%%%moving image%%%%%%%%%%%%%%%%" << endl;
> filter->GetOutput()->Print(std::cout);
> cout << endl;
>
> if (png) {
> typedef unsigned char OutputPixelType;
> typedef itk::Image<OutputPixelType, Dimension> OutputImageType;
> typedef itk::RescaleIntensityImageFilter<ImageType, OutputImageType>
> CastFilterType;
> CastFilterType::Pointer castFilter = CastFilterType::New();
> typedef itk::ImageFileWriter<OutputImageType> WriterType;
> WriterType::Pointer writer = WriterType::New();
> castFilter->SetOutputMinimum(0);
> castFilter->SetOutputMaximum(255);
> try {
> castFilter->SetInput(myImage);
> writer->SetFileName("d:\\regProjectData\\PNGs\\fixed.png");
> writer->SetInput(castFilter->GetOutput());
> writer->Update();
> castFilter->SetInput(filter->GetOutput());
> writer->SetFileName("d:\\regProjectData\\PNGs\\moving.png");
> writer->SetInput(castFilter->GetOutput());
> writer->Update();
> } catch (itk::ExceptionObject &err) {
> std::cout << "Exception while writing files: " << std::endl;
> std::cout << err << std::endl;
> return -2;
> }
> }
>
> registration->SetFixedImageRegion(
> myImage->GetBufferedRegion());
>
> typedef RegistrationType::ParametersType ParametersType;
> ParametersType initialParameters(
transform->GetNumberOfParameters());
> initialParameters[0] = FACTOR * 2;
> initialParameters[1] = FACTOR * 2;
> if (Dimension > 2) {
> initialParameters[2] = 0.0;
> }
> registration->SetInitialTransformParameters(initialParameters);
>
> double minStep;
> if (!vtkUtility::GetParam("-minStep", &minStep, argc, argv)) {
> minStep = FACTOR * 0.01;
> }
> cout << "using minimum step length of " << minStep << endl;
> double maxStep;
> if (!vtkUtility::GetParam("-maxStep", &maxStep, argc, argv)) {
> maxStep = FACTOR * 4.0;
> }
> cout << "using maximum step length of " << maxStep << endl;
> optimizer->SetMinimumStepLength(minStep);
> optimizer->SetMaximumStepLength(maxStep);
>
> optimizer->SetNumberOfIterations(200);
>
> Observer::Pointer observer = Observer::New();
>
> optimizer->AddObserver(itk::IterationEvent(), observer);
> optimizer->AddObserver(itk::StartEvent(), observer);
> optimizer->AddObserver(itk::EndEvent(), observer);
> optimizer->AddObserver(itk::ProgressEvent(), observer);
>
> try {
> registration->StartRegistration();
> } catch (itk::ExceptionObject &err) {
> std::cerr << "ExceptionObject caught!" << std::endl;
> std::cerr << err << std::endl;
> return -1;
> }
>
>
> ParametersType finalParameters =
> registration->GetLastTransformParameters();
> const double xTrans = finalParameters[0];
> const double yTrans = finalParameters[1];
> double zTrans = 0;
> if (Dimension > 2) {
> zTrans = finalParameters[2];
> }
>
> const unsigned int numberOfIterations =
optimizer->GetCurrentIteration();
> const double bestValue = optimizer->GetValue();
>
> std::cout << "xTrans: " << xTrans << std::endl;
> std::cout << "yTrans: " << yTrans << std::endl;
> if (Dimension > 2) {
> std::cout << "zTrans: " << zTrans << std::endl;
> }
>
> std::cout << "numberOfIterations: " << numberOfIterations <<
std::endl;
> std::cout << "bestValue: " << bestValue << std::endl;
>
>
> cout << "%%%%%%%%%%%%%%fixed image after reg%%%%%%%%%%%%%%%%%" << endl;
> myImage->Print(std::cout);
> cout << endl;
>
> //filter->Update();
>
> registration->SetMovingImage(filter->GetOutput());
> cout << "%%%%%%%%%%%%%%moving image after reg%%%%%%%%%%%%%%%%" <<
endl;
> filter->GetOutput()->Print(std::cout);
> cout << endl;
>
>
> return 0;
> }
>