[Insight-users] Re: SetFixedImageRegion
Wei-Lin Bee
tigerbee at MIT.EDU
Wed Sep 15 20:19:43 EDT 2004
Hi Luis,
Thank you very much for your instruction.
I've change MutualInformationImageToImageMetric to MeanSquaresImageToImageMetric
and replaced related parameters. Now the program seems work well.
One more question:
The output image seems to be resampled moving image at the time. How to combine
fixed image and moving image to a stitched one as output image?
Thanks abound!
Best,
Tiger
ps. again, the attachment is my program for your reference...
Quoting Luis Ibanez <luis.ibanez at kitware.com>:
>
> Hi Wei,
>
> Thanks for your detailed report.
>
> Your code seems to be in good shape.
>
> The run-time problem that you are facing is
> an exception being thrown when you invoke the
> method GetValue() in the MutualInformation metric.
>
> You *must* put GetValue() inside a try/catch block
> so you will be able to know when an exception is
> thrown. Please look at the ITK SoftwareGuide,
>
> http://www.itk.org/ItkSoftwareGuide.pdf
>
> there are multiple examples on how to catch exceptions.
>
>
> Once you add the try/catch block, this particular
> exception will show you the following message:
>
>
> itk::ExceptionObject (0x8199f10)
> Location: "Unknown"
> File:
>
/home/ibanez/src/Insight/Code/Algorithms/itkMutualInformationImageToImageMetric.txx
> Line: 155
> Description: itk::ERROR: MutualInformationImageToImageMetric(0x8139878):
> All the sampled point mapped to outside of the moving image
>
>
> Which means that MutualInformation was not able to find enough sampling
> point in the subregion in which you are performing this computation.
>
>
> As you probably know, after reading the Registration chapter
> from the SoftwareGuide, MutualInformation samples the image
> randomly in order to compute the joint histogram of gray levels.
> When not enough samples fall inside the FixedImageRegion, this
> exception is thrown, because the estimation is expected to be
> very poor with such a small number of samples.
>
>
> You could go around this by increasing the number of samples....
> but in practice, what we could suggest you is to NOT use Mutual
> Information as the metric for this task. This is based on two
> arguments:
>
>
> 1) Your images are of the same modality.
>
> 2) Your images are almost black with scattered bright signals
> (we assume they are cells in fluorescent microscopy).
> In that regard, information is only carried by the cells, not
> the background. The chances that some of the samples from Mutual
> Information will fall inside a cell are very small. In most cases
> your image will appear simply as a black field to the IM Metric.
>
>
> Please replace Mutual Information with a metric better suited for
> this problem, such as Normalized Correlation, or MeanSquares.
>
>
> BTW: Now you will not need to compute 100 times the metric and make
> an average from it. You will not need the Normalization filters
> either, since NormalizedCorrelation is invariant to linear changes
> in the intensity levels.
>
>
> Also, remember to put always *try/catch* blocks around the pieces
> of code that can potentially throw exceptions.
>
>
>
> Regards,
>
>
>
> Luis
>
>
>
> --------------------
> Wei-Lin Bee wrote:
> > Hi, Luis,
> >
> > These 10 days I hope I'd have done what you told me to do. Following
> > itkSoftwareGuide I made executable the attached program.
> >
> > However this program still has few questions:
> >
> > 1. RegionOfInterest problems
> >
> > I set region of interest 25% as follows:
> >
> > if(atoi(argv[4]) == 1)
> > {
> > startF[0] = 192;
> > startF[1] = 0;
> >
> > sizeF[0] = 64;
> > sizeF[1] = 256;
> >
> >
> > startM[0] = 0;
> > startM[1] = 0;
> >
> > sizeM[0] = 64;
> > sizeM[1] = 256;
> > }
> >
> >
> > else if(atoi(argv[4]) == 2)
> > {
> > startF[0] = 0;
> > startF[1] = 192;
> >
> > sizeF[0] = 256;
> > sizeF[1] = 64;
> >
> > startM[0] = 0;
> > startM[1] = 0;
> >
> > sizeM[0] = 256;
> > sizeM[1] = 64;
> > }
> >
> > else
> > {
> > startF[0] = 0;
> > startF[1] = 0;
> >
> > sizeF[0] = 256;
> > sizeF[1] = 256;
> >
> > startM[0] = 0;
> > startM[1] = 0;
> >
> > sizeM[0] = 256;
> > sizeM[1] = 256;
> > }
> >
> > --> When I chose 1 or 2, I got run time error:
> >
> > ex. two images,fixedImage and movingImage are in the same directory:
> >
> > C:>tigerbee fixed.mhd moving.mhd output.mhd 1
> >
> > This application has requested the Runtime to terminate it in an unusual
> way.
> > Please contact the application's support team for more information.
> >
> > --> The error appear when I chose start point != 0
> >
> > 2. Output Image Size:
> > After registration, I wanna get the whole stitched image. Everytime I
> got is
> > original size: 256 X 256 pixels, what I need should be 448 X 256 pixels.
> > I have no idea how to make it.
> >
> >
> >
> > Thank you very much.
> >
> >
> > Best,
> > Tiger
> > ps. Prevent misleading you, the attachments are my program and two image
> > files(moving image is to right of fixed image).
> >
> >
> >
> >
> > Quoting Luis Ibanez <luis.ibanez at kitware.com>:
> >
> >
> >>Hi Wei,
> >>
> >>
> >>Please do the following:
> >>
> >>
> >>1) Stop writing code.
> >>
> >>
> >>
> >>2) Read the section of resampling from the
> >> ITK Software Guide.
> >>
> >> http://www.itk.org/ItkSoftwareGuide.pdf
> >>
> >> Section 6.7.1, pdf-pages 199 to 218.
> >>
> >>
> >>
> >>3) Read the Chapter on Image Registration from
> >> the ITK SoftwareGuide. This is Chapter 8,
> >> pdf-pages 241 - 328.
> >>
> >> All the Examples available in
> >>
> >> Insight/Examples/Registration
> >>
> >> are explained in this Chapter of the software
> >> guide. To be more precise, the examples *are*
> >> the SoftwareGuide.
> >>
> >>
> >>
> >>Please let us know if you have any questions
> >>after reading this material,
> >>
> >>
> >> Thanks
> >>
> >>
> >> Luis
> >>
> >>
> >>
> >>--------------------
> >>Wei-Lin Bee wrote:
> >>
> >>
> >>>Hi Luis,
> >>>
> >>>You're right. My image files come from one modality. Sorry for misleading
> >>>you and thank you for correcting me immediately before I go too far.
> >>>
> >>>As for this case (roughly 5mm X 20mm overlap between 20mm X 20mm images),
>
> >>>Due to lack of generic programming knowledge, although you've already
> given
> >>
> >>me
> >>
> >>>clear direction and the next step I should follow, however, I don't
> >>
> >>exactly
> >>
> >>>know how to do.
> >>>
> >>>A.What I learned from API is:
> >>> In itk::ImageRegistrationMethod, there is SetFixedImageRegion function
> >>
> >>needed
> >>
> >>>to be instantiated?
> >>>
> >>>B.So I try to get similar code from Insight/Examples/Registration:
> >>>
> >>>1.ImageRegistration1~14:
> >>>
> >>> registration->SetFixedImageRegion(
> >>>fixedImageReader->GetOutput()->GetBufferedRegion() );
> >>>
> >>> or
> >>>
> >>> registration->SetFixedImageRegion(
> >>> fixedNormalizer->GetOutput()->GetBufferedRegion() );
> >>>
> >>>2.MeanSquaresImageMetric1
> >>>
> >>> metric->SetFixedImageRegion(fixedImage->GetBufferedRegion());
> >>>
> >>>3.MultiResImageRegistration1~2
> >>>
> >>> registration->SetFixedImageRegion(
> >>> fixedCaster->GetOutput()->GetBufferedRegion() );
> >>>
> >>>
> >>>4.DeformableRegistration6
> >>>
> >>> FixedImageType::RegionType fixedRegion =
> >>
> >>fixedImage->GetBufferedRegion();
> >>
> >>> registration->SetFixedImageRegion( fixedRegion );
> >>>
> >>>It seems that all of the examples are using full available content of
> >>
> >>images.
> >>
> >>>C. Finally, I seemingly got some hints from [Insight-users] mailing list
> >>>archives:
> >>>
> >>>
> >>
> >>ex.http://public.kitware.com/pipermail/insight-users/2004-March/007391.html
> >>
> >>> metric->SetMovingImageStandardDeviation(1.0);
> >>> metric->SetFixedImageStandardDeviation(1.0);
> >>> metric->SetFixedImageRegion(fixedImg->GetBufferedRegion());
> >>> metric->SetNumberOfSpatialSamples(1000);
> >>>
> >>>However, I think asking for a complete and accurate answer should be the
> >>
> >>best
> >>
> >>>way to learn ITK. (much better than aimless copy and paste myself)
> >>>On the other hand, I have wasted too much time(3 weeks)trial error in
> this
> >>
> >>test
> >>
> >>>sample registration(later, I have to apply ITK in 2D, 3D tissue images
> from
> >>
> >>16
> >>
> >>>channel two-photon microscopy). I got go faster.
> >>>
> >>>
> >>>
> >>>Thanks abound !
> >>>
> >>>
> >>>Best,
> >>>Tiger
> >>>
> >>>
> >>>
> >>>Quoting Luis Ibanez <luis.ibanez at kitware.com>:
> >>>
> >>>
> >>>
> >>>>Hi Tiger,
> >>>>
> >>>>
> >>>>1) It sounds like you are doing stitching in microscopy or
> >>>> creating a mosaic for video or satelite images.
> >>>>
> >>>> Is that the case ?
> >>>>
> >>>>
> >>>>2) The subject of your email says "Inter-Modality", but...
> >>>> if you are doning stitching, it is unlikely that your
> >>>> images are actually from multiple modalities...
> >>>>
> >>>> Please clarify:
> >>>>
> >>>> Are the images of the same modality or not ?
> >>>>
> >>>>
> >>>>3) If you are doing stiching you can solve the registration
> >>>> of every pair of images by setting the FixedImageRegion
> >>>> to the expected overlap (5mm x 20mm) and by making sure
> >>>> that you assign to all the images, values of pixel spacing
> >>>> and image origin that correspond roughly to their correct
> >>>> location in space. That will simplify a lot the registration.
> >>>>
> >>>>
> >>>>4) If you are doing this for Microscopy you may want to start
> >>>> with the following components:
> >>>>
> >>>> - NormalizedCorrelationMetric
> >>>> - TranslationTransform
> >>>> - LinearInterpolator
> >>>> - RegularStepGradientDescentOptimizer
> >>>>
> >>>>
> >>>>5) Make sure that you connect an Observer to the optimizer
> >>>> and you plot the values of the metric as the iterations
> >>>> proceed. That will guide in the process of fine tunning
> >>>> registration parameters. Keep in mind that you cannot
> >>>> trust an optimizer, you must supervise its progress and
> >>>> drive its parameters in order to obtain reasonable results.
> >>>> The use of command observers is illustrated in almost all
> >>>> the registration examples in
> >>>>
> >>>> Insight/Examples/Registration
> >>>>
> >>>>
> >>>>
> >>>>
> >>>> Regards,
> >>>>
> >>>>
> >>>> Luis
> >>>>
> >>>>
> >>>>
> >>>>-------------------------
> >>>>Wei-Lin Bee wrote:
> >>>>
> >>>>
> >>>>
> >>>>>Hi,
> >>>>>
> >>>>>I have 64 spacial image files (8 X 8)to be registered. The size of each
> >>>>
> >>>>image is
> >>>>
> >>>>
> >>>>>20mm X 20mm which overlap partion between two images would be 5mm X
> 20mm.
> >>>>
> >>>>I
> >>>>
> >>>>
> >>>>>learned from itkSoftwareGuide that there seems no clear-cut rules in
> >>>>
> >>>>choosing
> >>>>
> >>>>
> >>>>>matric, however I failed to solve this problem for days. Any reply
> would
> >>>>
> >>>>be
> >>>>
> >>>>
> >>>>>appreciated.
> >>>>>
> >>>>>
> >>>>>
> >>>>>Thank you very much.
> >>>>>
> >>>>>
> >>>>>Best,
> >>>>>Tiger
> >>>>>_______________________________________________
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>
> >>>
> >>
> >>
> >>
> >>
> >
> >
> >
> >
> > ------------------------------------------------------------------------
> >
> > //TigerBee.cxx
> >
> > #include "itkImageRegistrationMethod.h"
> > #include "itkTranslationTransform.h"
> > #include "itkMutualInformationImageToImageMetric.h"
> > #include "itkLinearInterpolateImageFunction.h"
> > #include "itkGradientDescentOptimizer.h"
> > #include "itkImage.h"
> >
> > #include "itkNormalizeImageFilter.h"
> > #include "itkImageFileReader.h"
> > #include "itkImageFileWriter.h"
> > #include "itkRegionOfInterestImageFilter.h"
> > #include "itkResampleImageFilter.h"
> > #include "itkCastImageFilter.h"
> > #include "itkSquaredDifferenceImageFilter.h"
> > #include <sstream>
> >
> > //
> > // The following section of code implements a Command observer
> > // that will monitor the evolution of the registration process.
> > //
> >
> > #include "itkCommand.h"
> > class CommandIterationUpdate : public itk::Command
> > {
> > public:
> > typedef CommandIterationUpdate Self;
> > typedef itk::Command Superclass;
> > typedef itk::SmartPointer<Self> Pointer;
> > itkNewMacro( Self );
> > protected:
> > CommandIterationUpdate() {};
> > public:
> > typedef itk::GradientDescentOptimizer OptimizerType;
> > typedef const OptimizerType * OptimizerPointer;
> >
> > void Execute(itk::Object *caller, const itk::EventObject & event)
> > {
> > Execute( (const itk::Object *)caller, event);
> > }
> >
> > void Execute(const itk::Object * object, const itk::EventObject & event)
> > {
> > #if 0
> > OptimizerPointer optimizer =
> > dynamic_cast< OptimizerPointer >( object );
> > if( typeid( event ) != typeid( itk::IterationEvent ) )
> > {
> > return;
> > }
> > std::cout << optimizer->GetCurrentIteration() << " ";
> > std::cout << optimizer->GetValue() << " ";
> > std::cout << optimizer->GetCurrentPosition()[0] << " ";
> > std::cout << optimizer->GetCurrentPosition()[1] << std::endl;
> > #endif
> > }
> > };
> >
> > int main( int argc, char **argv ){
> >
> > if( argc < 5 ){
> > std::cerr << "Missing Parameters " << std::endl;
> > std::cerr << "Usage: " << argv[0];
> > std::cerr << " fixedImageFile movingImageFile ";
> > std::cerr << " outputImagefile MovingImageLocation";
> > std::cerr << " [differenceOutputfile] [differenceBeforeRegistration]
> "<< std::endl;
> > return 1;
> > }
> >
> > // The types of images should be declared first.
> >
> > const unsigned int Dimension = 2;
> > typedef unsigned int PixelType;
> > typedef float InternalPixelType;
> > typedef float ROIPixelType;
> >
> > typedef itk::Image< PixelType, Dimension > FixedImageType;
> > typedef itk::Image< PixelType, Dimension > MovingImageType;
> > typedef itk::Image< ROIPixelType, Dimension > ROIFixedImageType;
> > typedef itk::Image< ROIPixelType, Dimension > ROIMovingImageType;
> > typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
> >
> > typedef itk::TranslationTransform< double, Dimension > TransformType;
> >
> > typedef itk::GradientDescentOptimizer OptimizerType;
> > typedef itk::LinearInterpolateImageFunction<
> > InternalImageType,
> > double > InterpolatorType;
> > typedef itk::MutualInformationImageToImageMetric<
> > InternalImageType,
> > InternalImageType >
> MetricType;
> >
> > TransformType::Pointer transform = TransformType::New();
> > InterpolatorType::Pointer interpolator = InterpolatorType::New();
> > MetricType::Pointer metric = MetricType::New();
> >
> > metric->SetFixedImageStandardDeviation( 0.4 );
> > metric->SetMovingImageStandardDeviation( 0.4 );
> > metric->SetNumberOfSpatialSamples( 100 );
> > metric->SetInterpolator( interpolator );
> > metric->SetTransform( transform );
> >
> > typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
> > typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
> >
> > // Settings for the region of interest filter for the moving and fixed
> images
> > typedef itk::RegionOfInterestImageFilter< FixedImageType,
> ROIFixedImageType > FixedROIFilterType;
> > typedef itk::RegionOfInterestImageFilter< MovingImageType,
> ROIMovingImageType > MovingROIFilterType;
> >
> > FixedROIFilterType::Pointer fixedROI = FixedROIFilterType::New();
> > MovingROIFilterType::Pointer movingROI = MovingROIFilterType::New();
> >
> > ROIFixedImageType::IndexType startF;
> > ROIFixedImageType::SizeType sizeF;
> > ROIMovingImageType::IndexType startM;
> > ROIMovingImageType::SizeType sizeM;
> >
> > // For MovingImageLocation (assume 25% overlap)
> > // 1 for moving image right of fixed image
> > // 2 for moving image above fixed image
> >
> > // If moving image is to right of fixed image
> >
> > if(atoi(argv[4]) == 1)
> > {
> > startF[0] = 192;
> > startF[1] = 0;
> >
> > sizeF[0] = 64;
> > sizeF[1] = 256;
> >
> >
> > startM[0] = 0;
> > startM[1] = 0;
> >
> > sizeM[0] = 64;
> > sizeM[1] = 256;
> > }
> >
> > // If moving image is above fixed image
> > else if(atoi(argv[4]) == 2)
> > {
> > startF[0] = 0;
> > startF[1] = 192;
> >
> > sizeF[0] = 256;
> > sizeF[1] = 64;
> >
> > startM[0] = 0;
> > startM[1] = 0;
> >
> > sizeM[0] = 256;
> > sizeM[1] = 64;
> > }
> > else
> > {
> > startF[0] = 0;
> > startF[1] = 0;
> >
> > sizeF[0] = 256;
> > sizeF[1] = 256;
> >
> > startM[0] = 0;
> > startM[1] = 0;
> >
> > sizeM[0] = 256;
> > sizeM[1] = 256;
> > }
> >
> > ROIFixedImageType::RegionType desiredFixedRegion;
> > desiredFixedRegion.SetIndex( startF );
> > desiredFixedRegion.SetSize( sizeF );
> >
> > ROIMovingImageType::RegionType desiredMovingRegion;
> > desiredMovingRegion.SetIndex( startM );
> > desiredMovingRegion.SetSize( sizeM );
> >
> > fixedROI->SetRegionOfInterest(desiredFixedRegion);
> > movingROI->SetRegionOfInterest(desiredMovingRegion);
> >
> > FixedImageReaderType::Pointer fixedImageReader =
> FixedImageReaderType::New();
> > MovingImageReaderType::Pointer movingImageReader =
> MovingImageReaderType::New();
> >
> > fixedImageReader->SetFileName( argv[1] );
> > movingImageReader->SetFileName( argv[2] );
> >
> > fixedROI->SetInput( fixedImageReader->GetOutput() );
> > movingROI->SetInput( movingImageReader->GetOutput() );
> >
> > // The normalization filters are declared using the fixed and moving
> image
> > // types as input and the internal image type as output.
> >
> > typedef itk::NormalizeImageFilter<
> > ROIFixedImageType, InternalImageType >
> FixedNormalizeFilterType;
> >
> > typedef itk::NormalizeImageFilter<
> > ROIMovingImageType, InternalImageType >
> MovingNormalizeFilterType;
> >
> > FixedNormalizeFilterType::Pointer fixedNormalizer =
> FixedNormalizeFilterType::New();
> >
> > MovingNormalizeFilterType::Pointer movingNormalizer =
> MovingNormalizeFilterType::New();
> >
> > // The output of the region of interest filters is connected as input to
> the
> > // normalization filters. The inputs to the registration method are
> taken
> > // from the normalization filters.
> >
> > fixedNormalizer->SetInput( fixedROI->GetOutput() );
> > movingNormalizer->SetInput( movingROI->GetOutput() );
> >
> > fixedNormalizer->Update();
> >
> > metric->SetFixedImage( fixedNormalizer->GetOutput() );
> > metric->SetMovingImage( movingNormalizer->GetOutput() );
> >
> > metric->SetFixedImageRegion(
> > fixedNormalizer->GetOutput()->GetBufferedRegion() );
> >
> > metric->Initialize();
> >
> > FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
> > FixedImageType::SizeType fixedSize =
> fixedImage->GetLargestPossibleRegion().GetSize();
> >
> > MetricType::TransformParametersType transformParameters( 2 );
> > transformParameters[0] = 0; //translation in x
> > transformParameters[1] = 0; //translation in y
> >
> > double best_translation_x = transformParameters[0];
> > double best_translation_y = transformParameters[1];
> > double best_mi_value = metric->GetValue( transformParameters );
> >
> > std::ofstream file_stream("Iterations.txt");
> >
> > for ( double x = 10.0; x <= 11.0; x += 0.125 ) {
> > for ( double y = -2; y <=0.0; y += 0.125 ) {
> > transformParameters[0] = x;
> > transformParameters[1] = y;
> > double sum = 0.0;
> > for ( int i = 0; i < 100; i++ ) {
> > double mi_value = metric->GetValue( transformParameters );
> > sum += mi_value;
> > #if 0
> > if ( mi_value > best_mi_value ) {
> > best_translation_x = x;
> > best_translation_y = y;
> > best_mi_value = mi_value;
> > }
> > #endif
> >
> > }
> > double avg = sum / 100.0;
> > std::cout << x << " " << y << " " << avg << '\n';
> > file_stream << x << " " << y << " " << avg << '\n';
> >
> > if ( avg > best_mi_value ) {
> > best_mi_value = avg;
> > best_translation_x = x;
> > best_translation_y = y;
> > }
> > }
> > }
> >
> >
> > double TranslationAlongX = best_translation_x;
> > double TranslationAlongY = best_translation_y;
> > double bestValue = best_mi_value;
> >
> > // Print out results
> > std::cout << "# Result = " << std::endl;
> > std::cout << "# Translation X = " << TranslationAlongX << std::endl;
> > std::cout << "# Translation Y = " << TranslationAlongY << std::endl;
> > std::cout << "# Metric value = " << bestValue << std::endl;
> >
> > file_stream << "# Translation X = " << TranslationAlongX <<'\n';
> > file_stream << "# Translation Y = " << TranslationAlongY <<'\n';
> > file_stream << "# Metric value = " << bestValue;
> > file_stream.close();
> >
> > typedef itk::ResampleImageFilter<
> > MovingImageType,
> > FixedImageType > ResampleFilterType;
> >
> > TransformType::Pointer finalTransform = TransformType::New();
> >
> > transformParameters[0] = best_translation_x+64; //translate transform
> plus amount of overlap
> > transformParameters[1] = best_translation_y;
> >
> > finalTransform->SetParameters( transformParameters );
> >
> > ResampleFilterType::Pointer resample = ResampleFilterType::New();
> >
> > resample->SetTransform( finalTransform );
> > resample->SetInput( movingImageReader->GetOutput() );
> > resample->SetSize( fixedSize );
> > /*
> > FixedImageType::SizeType size;
> > size[0]=448;
> > size[1]=256;
> > resample->SetSize(size);
> > */
> > resample->SetOutputOrigin( fixedImage->GetOrigin() );
> > resample->SetOutputSpacing( fixedImage->GetSpacing() );
> > resample->SetDefaultPixelValue( 100 );
> >
> >
> > typedef unsigned char OutputPixelType;
> >
> > typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
> >
> > typedef itk::CastImageFilter<
> > FixedImageType,
> > OutputImageType > CastFilterType;
> >
> > typedef itk::ImageFileWriter< OutputImageType > WriterType;
> >
> >
> > WriterType::Pointer writer = WriterType::New();
> > CastFilterType::Pointer caster = CastFilterType::New();
> >
> > writer->SetFileName( argv[3] );
> >
> > caster->SetInput( resample->GetOutput() );
> > writer->SetInput( caster->GetOutput() );
> > writer->Update();
> >
> > typedef itk::SquaredDifferenceImageFilter<
> > FixedImageType,
> > FixedImageType,
> > OutputImageType > DifferenceFilterType;
> >
> > DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
> >
> > WriterType::Pointer writer2 = WriterType::New();
> > writer2->SetInput( difference->GetOutput() );
> >
> > return 0;
> > }
>
>
>
>
>
-------------- next part --------------
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: TigerBee.cxx,v $
Language: C++
Date: $Date: 2004/09/15 20:21:36 $
Version: $Revision: 1.25 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#include "itkImageRegistrationMethod.h"
#include "itkTranslationTransform.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkRegionOfInterestImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"
int main( int argc, char *argv[] )
{
if( argc <5 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << " outputImagefile MovingImageLocation";
std::cerr << " [differenceOutputfile] [differenceBeforeRegistration] ";
return 1;
}
// The types of images should be declared first.
const unsigned int Dimension = 2;
typedef float PixelType;
typedef float InternalPixelType;
typedef float ROIPixelType;
typedef itk::Image< PixelType, Dimension > FixedImageType;
typedef itk::Image< PixelType, Dimension > MovingImageType;
typedef itk::Image< ROIPixelType, Dimension > ROIFixedImageType;
typedef itk::Image< ROIPixelType, Dimension > ROIMovingImageType;
typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
typedef itk::TranslationTransform< double, Dimension > TransformType;
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef itk::MeanSquaresImageToImageMetric< FixedImageType, MovingImageType > MetricType;
typedef itk:: LinearInterpolateImageFunction< MovingImageType, double > InterpolatorType;
typedef itk::ImageRegistrationMethod< FixedImageType, MovingImageType > RegistrationType;
MetricType::Pointer metric = MetricType::New();
TransformType::Pointer transform = TransformType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetMetric( metric );
registration->SetOptimizer( optimizer );
registration->SetTransform( transform );
registration->SetInterpolator( interpolator );
metric->SetInterpolator( interpolator );
metric->SetTransform( transform );
// Settings for the region of interest filter for the moving and fixed images
typedef itk::RegionOfInterestImageFilter<
FixedImageType,
ROIFixedImageType > FixedROIFilterType;
typedef itk::RegionOfInterestImageFilter<
MovingImageType,
ROIMovingImageType > MovingROIFilterType;
FixedROIFilterType::Pointer fixedROI = FixedROIFilterType::New();
MovingROIFilterType::Pointer movingROI = MovingROIFilterType::New();
ROIFixedImageType::IndexType startF;
ROIFixedImageType::SizeType sizeF;
ROIMovingImageType::IndexType startM;
ROIMovingImageType::SizeType sizeM;
// For MovingImageLocation (assume 40% overlap)
// 1 for moving image right of fixed image
// 2 for moving image above fixed image
// If moving image is to right of fixed image
if(atoi(argv[4]) == 1)
{
startF[0] = 154;
startF[1] = 0;
sizeF[0] = 102;
sizeF[1] = 256;
startM[0] = 0;
startM[1] = 0;
sizeM[0] = 154;
sizeM[1] = 256;
}
// If moving image is above fixed image
else if(atoi(argv[4]) == 2)
{
startF[0] = 0;
startF[1] = 154;
sizeF[0] = 256;
sizeF[1] = 102;
startM[0] = 0;
startM[1] = 0;
sizeM[0] = 256;
sizeM[1] = 102;
}
else
//if(atoi(argv[4]) == 3)
{
startF[0] = 0;
startF[1] = 0;
sizeF[0] = 256;
sizeF[1] = 256;
startM[0] = 0;
startM[1] = 0;
sizeM[0] = 256;
sizeM[1] = 256;
}
ROIFixedImageType::RegionType desiredFixedRegion;
desiredFixedRegion.SetIndex( startF );
desiredFixedRegion.SetSize( sizeF );
ROIMovingImageType::RegionType desiredMovingRegion;
desiredMovingRegion.SetIndex( startM );
desiredMovingRegion.SetSize( sizeM );
fixedROI->SetRegionOfInterest(desiredFixedRegion);
movingROI->SetRegionOfInterest(desiredMovingRegion);
FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
fixedImageReader->SetFileName( argv[1] );
movingImageReader->SetFileName( argv[2] );
fixedROI->SetInput( fixedImageReader->GetOutput() );
movingROI->SetInput( movingImageReader->GetOutput() );
try{
fixedROI->Update();
} catch (itk::ExceptionObject & err) {
std::cout<<"fixedROI ExceptionObject Caught!! "<<std::endl;
std::cout<<err<<std::endl;
}
metric->SetFixedImage( fixedImageReader->GetOutput() );
metric->SetMovingImage( movingImageReader->GetOutput() );
metric->SetFixedImageRegion(
fixedROI->GetOutput()->GetBufferedRegion() );
try{
metric->Initialize();
} catch ( itk::ExceptionObject & err ) {
std::cout<<"metric ExceptionObject Caught!!"<<std::endl;
std::cout<<err<<std::endl;
}
registration->SetFixedImage( fixedImageReader->GetOutput() );
registration->SetMovingImage( movingImageReader->GetOutput() );
try{
fixedImageReader->Update();
}catch(itk::ExceptionObject & err){
std::cout << " fixedImageReader ExceptionObject Caught!! " << std::endl;
std::cout << err <<std::endl;
}
registration->SetFixedImageRegion(
fixedROI->GetOutput()->GetBufferedRegion() );
typedef RegistrationType::ParametersType ParametersType;
ParametersType initialParameters( transform->GetNumberOfParameters() );
initialParameters[0] = 0.0; // Initial offset in mm along X
initialParameters[1] = 0.0; // Initial offset in mm along Y
registration->SetInitialTransformParameters( initialParameters );
optimizer->SetMaximumStepLength( 4.00 );
optimizer->SetMinimumStepLength( 0.01 );
optimizer->SetNumberOfIterations( 200 );
try
{
registration->Update();
} catch( itk::ExceptionObject & err ) {
std::cout << "registration ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return -1;
}
ParametersType finalParameters = registration->GetLastTransformParameters();
const double TranslationAlongX = finalParameters[0];
const double TranslationAlongY = finalParameters[1];
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const double bestValue = optimizer->GetValue();
std::cout << "Result = " << std::endl;
std::cout << " Translation X = " << TranslationAlongX << std::endl;
std::cout << " Translation Y = " << TranslationAlongY << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
typedef itk::ResampleImageFilter<
MovingImageType,
FixedImageType > ResampleFilterType;
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetInput( movingImageReader->GetOutput() );
resample->SetTransform( registration->GetOutput()->Get() );
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
// Tiger's question --> How to make stitched image?
//resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
FixedImageType::SizeType size;
size[0] = 500; // size along X
size[1] = 256; // size along Y
resample->SetSize(size);
//
resample->SetOutputOrigin( fixedImage->GetOrigin() );
resample->SetOutputSpacing( fixedImage->GetSpacing() );
resample->SetDefaultPixelValue( 100 );
typedef unsigned char OutputPixelType;
typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
typedef itk::CastImageFilter<
FixedImageType,
OutputImageType > CastFilterType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName( argv[3] );
caster->SetInput( resample->GetOutput() );
writer->SetInput( caster->GetOutput() );
try{
writer->Update();
}catch( itk::ExceptionObject & err ){
std::cout<<"writer ExceptionObject Caught!!"<<std::endl;
std::cout<< err <<std::endl;
}
typedef itk::SquaredDifferenceImageFilter<
FixedImageType,
FixedImageType,
OutputImageType > DifferenceFilterType;
DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
difference->SetInput1( fixedImageReader->GetOutput() );
difference->SetInput2( resample->GetOutput() );
WriterType::Pointer writer2 = WriterType::New();
writer2->SetInput( difference->GetOutput() );
if( argc >= 6 ){
difference->SetInput1( fixedImageReader->GetOutput() );
difference->SetInput2( resample->GetOutput() );
writer2->SetFileName( argv[5] );
try{
writer2->Update();
} catch ( itk::ExceptionObject & err ){
std::cout<<"argv[5] ExceptionObject Caught!!"<<std::endl;
std::cout<<err<<std::endl;
}
}
return 0;
}
More information about the Insight-users
mailing list