<DIV>Thanks Luis,</DIV>
<DIV>Sorry, but If you look again at the code you will see that I dit the initialization of CenteredTransformInitializer with the normalize and smooth images.</DIV>
<DIV>initializer->SetFixedImage( fixedSmoother->GetOutput() );<BR>initializer->SetMovingImage( movingSmoother->GetOutput() );</DIV>
<DIV>But if I replace these two lines with the orginals images:<BR>initializer->SetFixedImage( fixedImageReader->GetOutput() );<BR>initializer->SetMovingImage( movingImageReader->GetOutput() );<BR>I do not have an exception but it does not converge correctly.</DIV>
<DIV> </DIV>
<DIV> And I have to normalize and smooth my images to correctly use the Mutual information has you done in example 2....</DIV>
<DIV> </DIV>
<DIV>Just for information these are the value if I initialize without normalize and smooth the image. with the initial fixed and moving image.</DIV>
<DIV>Initial Transform Parameters<BR>[1, 0, 0, 1, 111.204, 131.591, -12.7911, -12.6914]<BR>0 0.360832 [0.993629, -0.00122832, 0.00831941, 1.01269, 111.204, 131.591, -12.7794, -12.5935]<BR>1 0.220862 [0.991423, -0.00269869, 0.00381388, 1.00791, 111.204, 131.591, -12.7428, -12.6268]<BR>2 0.458038 [0.991053, -0.00423641, 0.00555051, 1.0107, 111.204, 131.591, -12.7455, -12.6022]<BR>3 0.43253 [0.98915, -0.00600989, 0.00460683, 1.00976, 111.204, 131.591, -12.7564, -12.6076]<BR>4 0.391366 [0.988866, -0.0053906, 0.00362396, 1.0088, 111.204, 131.591, -12.7547, -12.6199]<BR>5 0.255813 [0.988056, -0.00601446, 0.00374833, 1.0088, 111.204, 131.591, -12.7609, -12.6203]<BR>6 0.23651 [0.987756, -0.00630565, 0.00309831, 1.00838, 111.204, 131.591, -12.7629, -12.6261]<BR>7 0.312396 [0.987953, -0.00558894,
0.00279363, 1.008, 111.204, 131.591, -12.7579, -12.6298]<BR>8 0.157922 [0.987289, -0.00548913, 0.00273997, 1.00797, 111.204, 131.591, -12.7529, -12.6261]<BR>9 0.259104 [0.986384, -0.0058869, 0.00200088, 1.00744, 111.204, 131.591, -12.7518, -12.6287]<BR>10 0.37124 [0.98628, -0.00594178, 0.00209053, 1.00749, 111.204, 131.591, -12.753, -12.6278]<BR>Result =<BR> Center X = 111.204<BR> Center Y = 131.591<BR> Translation X = -12.753<BR> Translation Y = -12.6278<BR> Iterations = 12<BR> Metric value = 0.520743</DIV>
<DIV> </DIV>
<DIV> </DIV>
<DIV> </DIV>
<DIV>thanks again for your help</DIV>
<DIV>Lagaffe<BR><BR><B><I>Luis Ibanez <luis.ibanez@kitware.com></I></B> wrote:</DIV>
<BLOCKQUOTE class=replbq style="PADDING-LEFT: 5px; MARGIN-LEFT: 5px; BORDER-LEFT: #1010ff 2px solid"><BR><BR>Hi Lagaffe,<BR><BR><BR>The error is likely related to the fact that you commented<BR>out the lines where the input images are connected to the <BR>CenteredTransformInitializer<BR><BR>><BR>> /*<BR>> initializer->SetFixedImage( fixedImageReader->GetOutput() );<BR>> initializer->SetMovingImage( movingImageReader->GetOutput() );<BR>> */<BR>><BR><BR><BR>Therefore, the initialization of the Transform is not being<BR>performed correctly and you start your optimization process<BR>with corrupted numerical values such as:<BR><BR><BR>> ---------------<BR>> Initial Transform Parameters<BR>> [1, 0, 0, 1, -1.00276e+008, -4.25387e+008, 4.90708e+007, 3.85172e+008]<BR>><BR><BR><BR>Also, note that you *MUST* call Update() in the readers before you use<BR>the TransformInitializer object. The reason is that the Initializer is<BR>not a pipeline object and
therefore it will not trigger an update in its<BR>inputs. This object is similar to the image calculator in this regard.<BR><BR><BR>If you notice that the optimizer is stopping too early without<BR>converging, that indicates that you are using:<BR><BR>A) a value too large for the step length of the optimizer<BR>B) inappropriate values for the parameter scaling<BR><BR><BR>If you continue having problems while fine tunning those parameters<BR>please post to the list the log output of the Command/Observer.<BR>From the plot of the Metric values and the Transform parameters at<BR>every iteration it is usually easier to figure out what parameters<BR>must be changed in your optimization setup.<BR><BR><BR><BR><BR>Please let us know if you have further questions,<BR><BR><BR><BR>Thanks,<BR><BR><BR><BR>Luis<BR><BR><BR><BR>---------------<BR>Mr Gaffe wrote:<BR><BR>> Hi all,<BR>> I try to combine registration example: ImageRegistration2 and 9 to <BR>> perform an
CenteredAffine+MutualInfo registration.<BR>> The code bellow compile well, but when I try with: <BR>> BrainT1SliceBorder20.png and BrainProtonDensitySliceR10X13Y17S12.png I <BR>> have the following error:<BR>> <BR>> ---------------<BR>> Initial Transform Parameters<BR>> [1, 0, 0, 1, -1.00276e+008, -4.25387e+008, 4.90708e+007, 3.85172e+008]<BR>> ExceptionObject caught !<BR>> itk::ExceptionObject (00FCFD30)<BR>> Location: "Unknown"<BR>> File: <BR>> D:\InsightToolkit-1.8.1\Code\Algorithms\itkMutualInformationImageToImageMetric.txx<BR>> Line: 155<BR>> Description: itk::ERROR: MutualInformationImageToImageMetric(012F4780): <BR>> All the sampled point mapped to outside of the moving image<BR>> -----------------<BR>> <BR>> ->The Typedef is set to float since Images are unsigned chart is it the <BR>> problem ?<BR>> If I change the initialisation of the CenteredTransformInitializer as:<BR>> initializer->SetFixedImage(
fixedImageReader->GetOutput() );<BR>> initializer->SetMovingImage( movingImageReader->GetOutput() );<BR>> There is no error, but the after 9 iterations it stop with a wrong image <BR>> results....<BR>> <BR>> Is it just a problem of initial typeddef for my fixed and moved image ?<BR>> <BR>> thanks,<BR>> <BR>> lagaffe<BR>> <BR>> ---------------------------------------------------------------------------------------------------------------------<BR>> /*<BR>> * Code from ImageRegistration2 and ImageRegistration9<BR>> */<BR>> #include "itkImageRegistrationMethod.h"<BR>> //#include "itkMeanSquaresImageToImageMetric.h"<BR>> #include "itkMutualInformationImageToImageMetric.h"<BR>> #include "itkNormalizeImageFilter.h"<BR>> #include "itkDiscreteGaussianImageFilter.h"<BR>> #include "itkLinearInterpolateImageFunction.h"<BR>> #include "itkRegularStepGradientDescentOptimizer.h"<BR>> #include "itkImage.h"<BR>> <BR>>
#include "itkCenteredTransformInitializer.h"<BR>> #include "itkCenteredAffineTransform.h"<BR>> #include "itkImageFileReader.h"<BR>> #include "itkImageFileWriter.h"<BR>> #include "itkResampleImageFilter.h"<BR>> #include "itkCastImageFilter.h"<BR>> #include "itkSquaredDifferenceImageFilter.h"<BR>> #include "itkCommand.h"<BR>> class CommandIterationUpdate : public itk::Command<BR>> {<BR>> public:<BR>> typedef CommandIterationUpdate Self;<BR>> typedef itk::Command Superclass;<BR>> typedef itk::SmartPointer<SELF> Pointer;<BR>> itkNewMacro( Self );<BR>> protected:<BR>> CommandIterationUpdate() {};<BR>> public:<BR>> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;<BR>> typedef const OptimizerType * OptimizerPointer;<BR>> void Execute(itk::Object *caller, const itk::EventObject & event)<BR>> {<BR>> Execute( (const itk::Object *)caller, event);<BR>> }<BR>> void Execute(const itk::Object * object, const
itk::EventObject & event)<BR>> {<BR>> OptimizerPointer optimizer =<BR>> dynamic_cast< OptimizerPointer >( object );<BR>> if( typeid( event ) != typeid( itk::IterationEvent ) )<BR>> {<BR>> return;<BR>> }<BR>> std::cout << optimizer->GetCurrentIteration() << " ";<BR>> std::cout << optimizer->GetValue() << " ";<BR>> std::cout << optimizer->GetCurrentPosition() << std::endl;<BR>> }<BR>> };<BR>> <BR>> int main( int argc, char *argv[] )<BR>> {<BR>> if( argc < 4 )<BR>> {<BR>> std::cerr << "Missing Parameters " << std::endl;<BR>> std::cerr << "Usage: " << argv[0];<BR>> std::cerr << " fixedImageFile movingImageFile " << std::endl;<BR>> std::cerr << " outputImagefile [differenceOutputfile] <BR>> [differenceBeforeRegistration] " << std::endl;<BR>> std::cerr << " [stepLength] [maxNumberOfIterations] "<<
std::endl;<BR>> return 1;<BR>> }<BR>> const unsigned int Dimension = 2;<BR>> typedef float PixelType;<BR>> typedef float InternalPixelType;<BR>> typedef itk::Image< InternalPixelType, Dimension > InternalImageType;<BR>> typedef itk::Image< PixelType, Dimension > FixedImageType;<BR>> typedef itk::Image< PixelType, Dimension > MovingImageType;<BR>> <BR>> typedef itk::CenteredAffineTransform<<BR>> double,<BR>> Dimension > TransformType;<BR>> typedef itk::RegularStepGradientDescentOptimizer OptimizerType;<BR>> /*<BR>> typedef itk::MeanSquaresImageToImageMetric<<BR>> FixedImageType,<BR>> MovingImageType > MetricType;<BR>> */<BR>> typedef itk::MutualInformationImageToImageMetric<<BR>> InternalImageType,<BR>> InternalImageType > MetricType;<BR>> typedef itk:: LinearInterpolateImageFunction<<BR>> MovingImageType,<BR>> double > InterpolatorType;<BR>> typedef
itk::ImageRegistrationMethod<<BR>> FixedImageType,<BR>> MovingImageType > RegistrationType;<BR>> MetricType::Pointer metric = MetricType::New();<BR>> OptimizerType::Pointer optimizer = OptimizerType::New();<BR>> InterpolatorType::Pointer interpolator = InterpolatorType::New();<BR>> RegistrationType::Pointer registration = RegistrationType::New();<BR>> TransformType::Pointer transform = TransformType::New();<BR>> <BR>> registration->SetMetric( metric );<BR>> registration->SetOptimizer( optimizer );<BR>> registration->SetInterpolator( interpolator );<BR>> registration->SetTransform( transform );<BR>> metric->SetFixedImageStandardDeviation( 0.4 );<BR>> metric->SetMovingImageStandardDeviation( 0.4 );<BR>> metric->SetNumberOfSpatialSamples( 50 );<BR>> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;<BR>> typedef itk::ImageFileReader< MovingImageType >
MovingImageReaderType;<BR>> FixedImageReaderType::Pointer fixedImageReader = <BR>> FixedImageReaderType::New();<BR>> MovingImageReaderType::Pointer movingImageReader = <BR>> MovingImageReaderType::New();<BR>> fixedImageReader->SetFileName( argv[1] );<BR>> movingImageReader->SetFileName( argv[2] );<BR>> typedef itk::NormalizeImageFilter<<BR>> FixedImageType,<BR>> InternalImageType<BR>> > FixedNormalizeFilterType;<BR>> typedef itk::NormalizeImageFilter<<BR>> MovingImageType,<BR>> InternalImageType<BR>> > MovingNormalizeFilterType;<BR>> FixedNormalizeFilterType::Pointer fixedNormalizer =<BR>> FixedNormalizeFilterType::New();<BR>> MovingNormalizeFilterType::Pointer movingNormalizer =<BR>> <BR>> MovingNormalizeFilterType::New();<BR>> <BR>> typedef itk::DiscreteGaussianImageFilter<<BR>> InternalImageType,<BR>> InternalImageType<BR>> > GaussianFilterType;<BR>> <BR>>
GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();<BR>> GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();<BR>> fixedSmoother->SetVariance( 2.0 );<BR>> movingSmoother->SetVariance( 2.0 );<BR>> <BR>> fixedNormalizer->SetInput( fixedImageReader->GetOutput() );<BR>> movingNormalizer->SetInput( movingImageReader->GetOutput() );<BR>> fixedSmoother->SetInput( fixedNormalizer->GetOutput() );<BR>> movingSmoother->SetInput( movingNormalizer->GetOutput() );<BR>> registration->SetFixedImage( fixedSmoother->GetOutput() );<BR>> registration->SetMovingImage( movingSmoother->GetOutput() );<BR>> <BR>> fixedNormalizer->Update();<BR>> registration->SetFixedImageRegion(<BR>> fixedNormalizer->GetOutput()->GetBufferedRegion() );<BR>> <BR>> typedef itk::CenteredTransformInitializer<<BR>> TransformType,<BR>> FixedImageType,<BR>> MovingImageType >
<BR>> TransformInitializerType;<BR>> TransformInitializerType::Pointer initializer = <BR>> TransformInitializerType::New();<BR>> initializer->SetTransform( transform );<BR>> /*<BR>> initializer->SetFixedImage( fixedImageReader->GetOutput() );<BR>> initializer->SetMovingImage( movingImageReader->GetOutput() );<BR>> */<BR>> <BR>> initializer->SetFixedImage( fixedSmoother->GetOutput() );<BR>> initializer->SetMovingImage( movingSmoother->GetOutput() );<BR>> initializer->MomentsOn();<BR>> initializer->InitializeTransform();<BR>> std::cout << "Initial Transform Parameters " << std::endl;<BR>> std::cout << transform->GetParameters() << std::endl;<BR>> <BR>> // -----------------------<BR>> registration->SetInitialTransformParameters(<BR>> transform->GetParameters() );<BR>> <BR>> double translationScale = 1.0 / 1000.0;<BR>> if( argc > 8 )<BR>> {<BR>>
translationScale = atof( argv[8] );<BR>> }<BR>> typedef OptimizerType::ScalesType OptimizerScalesType;<BR>> OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );<BR>> optimizerScales[0] = 1.0;<BR>> optimizerScales[1] = 1.0;<BR>> optimizerScales[2] = 1.0;<BR>> optimizerScales[3] = 1.0;<BR>> optimizerScales[4] = translationScale;<BR>> optimizerScales[5] = translationScale;<BR>> optimizerScales[6] = translationScale;<BR>> optimizerScales[7] = translationScale;<BR>> optimizer->SetScales( optimizerScales );<BR>> double steplength = 0.1;<BR>> if( argc > 6 )<BR>> {<BR>> steplength = atof( argv[6] );<BR>> }<BR>> <BR>> unsigned int maxNumberOfIterations = 300;<BR>> if( argc > 7 )<BR>> {<BR>> maxNumberOfIterations = atoi( argv[7] );<BR>> }<BR>> optimizer->SetMaximumStepLength( steplength );<BR>> optimizer->SetMinimumStepLength( 0.001 );<BR>> optimizer->SetNumberOfIterations(
maxNumberOfIterations );<BR>> optimizer->MinimizeOn();<BR>> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();<BR>> optimizer->AddObserver( itk::IterationEvent(), observer );<BR>> try<BR>> {<BR>> registration->StartRegistration();<BR>> }<BR>> catch( itk::ExceptionObject & err )<BR>> {<BR>> std::cerr << "ExceptionObject caught !" << std::endl;<BR>> std::cerr << err << std::endl;<BR>> return -1;<BR>> }<BR>> OptimizerType::ParametersType finalParameters =<BR>> registration->GetLastTransformParameters();<BR>> const double finalRotationCenterX = finalParameters[4];<BR>> const double finalRotationCenterY = finalParameters[5];<BR>> const double finalTranslationX = finalParameters[6];<BR>> const double finalTranslationY = finalParameters[7];<BR>> const unsigned int numberOfIterations = optimizer->GetCurrentIteration();<BR>> const double bestValue =
optimizer->GetValue();<BR>> // Software Guide : EndCodeSnippet<BR>> <BR>> // Print out results<BR>> //<BR>> std::cout << "Result = " << std::endl;<BR>> std::cout << " Center X = " << finalRotationCenterX << std::endl;<BR>> std::cout << " Center Y = " << finalRotationCenterY << std::endl;<BR>> std::cout << " Translation X = " << finalTranslationX << std::endl;<BR>> std::cout << " Translation Y = " << finalTranslationY << std::endl;<BR>> std::cout << " Iterations = " << numberOfIterations << std::endl;<BR>> std::cout << " Metric value = " << bestValue << std::endl;<BR>> typedef itk::ResampleImageFilter<<BR>> MovingImageType,<BR>> FixedImageType > ResampleFilterType;<BR>> TransformType::Pointer finalTransform = TransformType::New();<BR>> finalTransform->SetParameters( finalParameters );<BR>>
ResampleFilterType::Pointer resample = ResampleFilterType::New();<BR>> resample->SetTransform( finalTransform );<BR>> resample->SetInput( movingImageReader->GetOutput() );<BR>> FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();<BR>> resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );<BR>> resample->SetOutputOrigin( fixedImage->GetOrigin() );<BR>> resample->SetOutputSpacing( fixedImage->GetSpacing() );<BR>> resample->SetDefaultPixelValue( 100 );<BR>> <BR>> typedef unsigned char OutputPixelType;<BR>> typedef itk::Image< OutputPixelType, Dimension > OutputImageType;<BR>> <BR>> typedef itk::CastImageFilter<<BR>> FixedImageType,<BR>> OutputImageType > CastFilterType;<BR>> <BR>> typedef itk::ImageFileWriter< OutputImageType > WriterType;<BR>> <BR>> WriterType::Pointer writer = WriterType::New();<BR>> CastFilterType::Pointer caster =
CastFilterType::New();<BR>> <BR>> writer->SetFileName( argv[3] );<BR>> <BR>> caster->SetInput( resample->GetOutput() );<BR>> writer->SetInput( caster->GetOutput() );<BR>> writer->Update();<BR>> <BR>> typedef itk::SquaredDifferenceImageFilter<<BR>> FixedImageType,<BR>> FixedImageType,<BR>> OutputImageType > DifferenceFilterType;<BR>> DifferenceFilterType::Pointer difference = DifferenceFilterType::New();<BR>> WriterType::Pointer writer2 = WriterType::New();<BR>> writer2->SetInput( difference->GetOutput() ); <BR>> <BR>> if( argc > 4 )<BR>> {<BR>> difference->SetInput1( fixedImageReader->GetOutput() );<BR>> difference->SetInput2( resample->GetOutput() );<BR>> writer2->SetFileName( argv[4] );<BR>> writer2->Update();<BR>> }<BR>> <BR>> if( argc > 5 )<BR>> {<BR>> writer2->SetFileName( argv[5] );<BR>> difference->SetInput1(
fixedImageReader->GetOutput() );<BR>> difference->SetInput2( movingImageReader->GetOutput() );<BR>> writer2->Update();<BR>> }<BR>> <BR>> return 0;<BR>> }<BR>> <BR>> <BR>> ------------------------------------------------------------------------<BR>> Créez gratuitement votre Yahoo! Mail avec *100 Mo de stockage !*<BR>> Créez votre Yahoo! Mail <BR>> <HTTP: evt="25917/*http://fr.rd.yahoo.com/mail/mail_taglines_100/default/*http://fr.benefits.yahoo.com/" fr.rd.yahoo.com *http: taglines mail><BR>> <BR>> *Le nouveau Yahoo! Messenger est arrivé !* Découvrez toutes les <BR>> nouveautés pour dialoguer instantanément avec vos amis. Téléchargez <BR>> GRATUITEMENT ici ! <BR>> <HTTP: evt="26111/*http://fr.rd.yahoo.com/messenger/mail_taglines/default/*http://fr.messenger.yahoo.com" fr.rd.yahoo.com *http: taglines mail><BR>> <BR>> <BR>> <BR>> ------------------------------------------------------------------------<BR>>
<BR>> _______________________________________________<BR>> Insight-users mailing list<BR>> Insight-users@itk.org<BR>> http://www.itk.org/mailman/listinfo/insight-users<BR><BR><BR><BR><BR></BLOCKQUOTE><p>
<hr size=1>
Créez gratuitement votre Yahoo! Mail avec <font color="red"><b>100 Mo de stockage !</b></font>
<br><a href="http://fr.rd.yahoo.com/mail/taglines/*http://fr.rd.yahoo.com/evt=25917/*http://fr.rd.yahoo.com/mail/mail_taglines_100/default/*http://fr.benefits.yahoo.com/">Créez votre Yahoo! Mail</a><br><br>
<font color="red"><b>Le nouveau Yahoo! Messenger est arrivé !</b></font> Découvrez toutes les nouveautés pour dialoguer instantanément avec vos amis.
<a href="http://fr.rd.yahoo.com/mail/taglines/*http://fr.rd.yahoo.com/evt=26111/*http://fr.rd.yahoo.com/messenger/mail_taglines/default/*http://fr.messenger.yahoo.com">Téléchargez GRATUITEMENT ici !</a>