[Insight-users] Re: SetFixedImageRegion
Luis Ibanez
luis.ibanez at kitware.com
Sat Sep 11 13:30:09 EDT 2004
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;
> }
More information about the Insight-users
mailing list