Hi All:<br><br>After trial and error in my program that include Mattes MI, Affine transform, and RegularStepGradientDescentOptimizer.<br>The result is also appearing dispiritedly.<br><br>I have several questions and appreciated with your help.
<br><br>1)Is it right to simply set Mattes MI using <br>&nbsp; metric-&gt;SetNumberOfHistogramBins( 50 );<br>&nbsp; metric-&gt;SetNumberOfSpatialSamples( 10000 );<br>&nbsp; metric-&gt;ReinitializeSeed( 76926294 );&nbsp; ?<br><br>2)Is it right to set Optimizer with 
<br>&nbsp; optimizerScales[0] =&nbsp; 1.0;<br>&nbsp; optimizerScales[1] =&nbsp; 1.0;<br>&nbsp; optimizerScales[2] =&nbsp; 1.0;<br>&nbsp; optimizerScales[3] =&nbsp; 1.0;<br>&nbsp; optimizerScales[4] =&nbsp; 1/1000000;<br>&nbsp; optimizerScales[5] =&nbsp; 1/1000000;<br>&nbsp; optimizer-&gt;SetScales( optimizerScales );
<br>&nbsp; optimizer-&gt;MinimizeOn();<br>&nbsp; optimizer-&gt;SetMaximumStepLength( 0.5 ); <br>&nbsp; optimizer-&gt;SetMinimumStepLength( 0.1 );<br>&nbsp; optimizer-&gt;SetNumberOfIterations( 200 );<br>&nbsp; optimizer-&gt;SetRelaxationFactor( 0.8
 );<br><br>3) How to determine an appropriate parameter for Optimizer?<br><br>My program is following as below, If anyone can give me any suggestion I will appreciate with your help.<br><br>Goo :)<br><br>=======================================================================
<br><br>int main( int argc, char *argv[] )<br>{<br>&nbsp; if( argc &lt; 4 )<br>&nbsp;&nbsp;&nbsp; {<br>&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; &quot;Missing Parameters &quot; &lt;&lt; std::endl;<br>&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; &quot;Usage: &quot; &lt;&lt; argv[0];
<br>&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; &quot;&nbsp;&nbsp; fixedImageFile&nbsp; movingImageFile &quot; &lt;&lt; std::endl;<br>&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; &quot;&nbsp;&nbsp; outputImagefile&nbsp; [differenceBeforeRegistration] &quot; &lt;&lt; std::endl;<br>&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; &quot;&nbsp;&nbsp; [differenceAfterRegistration] &quot; &lt;&lt; std::endl;
<br>&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; &quot;&nbsp;&nbsp; [stepLength] [maxNumberOfIterations] &quot;&lt;&lt; std::endl;<br>&nbsp;&nbsp;&nbsp; return 1;<br>&nbsp;&nbsp;&nbsp; }<br><br>&nbsp; const&nbsp;&nbsp;&nbsp; unsigned int&nbsp;&nbsp;&nbsp; Dimension = 2;<br>&nbsp; typedef&nbsp; float&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; PixelType;<br><br>
&nbsp; typedef itk::Image&lt; PixelType, Dimension &gt;&nbsp; FixedImageType;<br>&nbsp; typedef itk::Image&lt; PixelType, Dimension &gt;&nbsp; MovingImageType;<br><br>&nbsp; typedef itk::AffineTransform&lt; double, Dimension &gt;&nbsp;&nbsp;&nbsp;&nbsp; TransformType;
<br><br>&nbsp; typedef itk::RegularStepGradientDescentOptimizer&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OptimizerType;<br><br><br>&nbsp; typedef itk::MattesMutualInformationImageToImageMetric&lt; <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; FixedImageType, <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; MovingImageType &gt;&nbsp;&nbsp;&nbsp; MetricType;
<br><br><br>&nbsp; typedef itk:: LinearInterpolateImageFunction&lt; <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; MovingImageType,<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; double&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &gt;&nbsp;&nbsp;&nbsp; InterpolatorType;<br>&nbsp; <br>&nbsp; typedef itk::ImageRegistrationMethod&lt; 
<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; FixedImageType, <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; MovingImageType &gt;&nbsp;&nbsp;&nbsp; RegistrationType;<br><br>&nbsp; MetricType::Pointer&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; metric&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = MetricType::New();<br>&nbsp; OptimizerType::Pointer&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; optimizer&nbsp;&nbsp;&nbsp;&nbsp; = OptimizerType::New();
<br>&nbsp; InterpolatorType::Pointer&nbsp;&nbsp; interpolator&nbsp; = InterpolatorType::New();<br>&nbsp; RegistrationType::Pointer&nbsp;&nbsp; registration&nbsp; = RegistrationType::New();<br><br>&nbsp; registration-&gt;SetMetric(&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; metric&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; );<br>&nbsp; registration-&gt;SetOptimizer(&nbsp;&nbsp;&nbsp;&nbsp; optimizer&nbsp;&nbsp;&nbsp;&nbsp; );
<br>&nbsp; registration-&gt;SetInterpolator(&nbsp; interpolator&nbsp; );<br><br>&nbsp; TransformType::Pointer&nbsp; transform = TransformType::New();<br>&nbsp; registration-&gt;SetTransform( transform );<br><br><br>&nbsp; metric-&gt;SetNumberOfHistogramBins( 50 );
<br>&nbsp; metric-&gt;SetNumberOfSpatialSamples( 10000 );<br>&nbsp; metric-&gt;ReinitializeSeed( 76926294 );<br><br><br>&nbsp; typedef itk::ImageFileReader&lt; FixedImageType&nbsp; &gt; FixedImageReaderType;<br>&nbsp; typedef itk::ImageFileReader&lt; MovingImageType &gt; MovingImageReaderType;
<br>&nbsp; FixedImageReaderType::Pointer&nbsp; fixedImageReader&nbsp; = FixedImageReaderType::New();<br>&nbsp; MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();<br>&nbsp; fixedImageReader-&gt;SetFileName(&nbsp; argv[1] );
<br>&nbsp; movingImageReader-&gt;SetFileName( argv[2] );<br><br><br>&nbsp; registration-&gt;SetFixedImage(&nbsp;&nbsp;&nbsp; fixedImageReader-&gt;GetOutput()&nbsp;&nbsp;&nbsp; );<br>&nbsp; registration-&gt;SetMovingImage(&nbsp;&nbsp; movingImageReader-&gt;GetOutput()&nbsp;&nbsp; );<br>
<br>&nbsp; fixedImageReader-&gt;Update();<br><br>&nbsp; registration-&gt;SetFixedImageRegion( <br>&nbsp;&nbsp;&nbsp;&nbsp; fixedImageReader-&gt;GetOutput()-&gt;GetBufferedRegion() );<br><br><br>&nbsp; typedef itk::CenteredTransformInitializer&lt; <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; TransformType, 
<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; FixedImageType, <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; MovingImageType &gt;&nbsp; TransformInitializerType;<br>&nbsp; TransformInitializerType::Pointer initializer = TransformInitializerType::New();
<br>&nbsp; initializer-&gt;SetTransform(&nbsp;&nbsp; transform );<br>&nbsp; initializer-&gt;SetFixedImage(&nbsp; fixedImageReader-&gt;GetOutput() );<br>&nbsp; initializer-&gt;SetMovingImage( movingImageReader-&gt;GetOutput() );<br>&nbsp; initializer-&gt;MomentsOn();
<br>&nbsp; initializer-&gt;InitializeTransform();<br><br>&nbsp; registration-&gt;SetInitialTransformParameters( <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; transform-&gt;GetParameters() );<br><br><br>&nbsp; double translationScale = 1.0 / 100000.0
;<br>&nbsp; if( argc &gt; 8 ){<br>&nbsp;&nbsp;&nbsp; translationScale = atof( argv[8] );<br>&nbsp;&nbsp;&nbsp; }<br><br>&nbsp; typedef OptimizerType::ScalesType&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OptimizerScalesType;<br>&nbsp; OptimizerScalesType optimizerScales( transform-&gt;GetNumberOfParameters() );
<br><br>&nbsp; optimizerScales[0] =&nbsp; 1.0;<br>&nbsp; optimizerScales[1] =&nbsp; 1.0;<br>&nbsp; optimizerScales[2] =&nbsp; 1.0;<br>&nbsp; optimizerScales[3] =&nbsp; 1.0;<br>&nbsp; optimizerScales[4] =&nbsp; 1/1000000;<br>&nbsp; optimizerScales[5] =&nbsp; 1/1000000;<br><br>&nbsp; optimizer-&gt;SetScales( optimizerScales );
<br><br><br>&nbsp; optimizer-&gt;MinimizeOn();<br>&nbsp; optimizer-&gt;SetMaximumStepLength( 0.5 ); <br>&nbsp; optimizer-&gt;SetMinimumStepLength( 0.1 );<br>&nbsp; optimizer-&gt;SetNumberOfIterations( 200 );<br>&nbsp; optimizer-&gt;SetRelaxationFactor( 
0.8 );<br>&nbsp; <br><br>&nbsp; CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();<br>&nbsp; optimizer-&gt;AddObserver( itk::IterationEvent(), observer );<br><br>&nbsp; try <br>&nbsp;&nbsp;&nbsp; { <br>&nbsp;&nbsp;&nbsp; registration-&gt;StartRegistration(); 
<br>&nbsp;&nbsp;&nbsp; } <br>&nbsp; catch( itk::ExceptionObject &amp; err ) <br>&nbsp;&nbsp;&nbsp; { <br>&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; &quot;ExceptionObject caught !&quot; &lt;&lt; std::endl; <br>&nbsp;&nbsp;&nbsp; std::cerr &lt;&lt; err &lt;&lt; std::endl; <br>&nbsp;&nbsp;&nbsp; return -1;<br>
&nbsp;&nbsp;&nbsp; } <br><br>&nbsp; OptimizerType::ParametersType finalParameters = <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; registration-&gt;GetLastTransformParameters();<br><br>&nbsp; const double finalRotationCenterX = transform-&gt;GetCenter()[0];<br>&nbsp; const double finalRotationCenterY = transform-&gt;GetCenter()[1];
<br>&nbsp; const double finalTranslationX&nbsp;&nbsp;&nbsp; = finalParameters[4];<br>&nbsp; const double finalTranslationY&nbsp;&nbsp;&nbsp; = finalParameters[5];<br><br>&nbsp; const unsigned int numberOfIterations = optimizer-&gt;GetCurrentIteration();<br>&nbsp; const double bestValue = optimizer-&gt;GetValue();
<br>&nbsp;<br>&nbsp; std::cout &lt;&lt; &quot;Result = &quot; &lt;&lt; std::endl;<br>&nbsp; std::cout &lt;&lt; &quot; Center X&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = &quot; &lt;&lt; finalRotationCenterX&nbsp; &lt;&lt; std::endl;<br>&nbsp; std::cout &lt;&lt; &quot; Center Y&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = &quot; &lt;&lt; finalRotationCenterY&nbsp; &lt;&lt; std::endl;
<br>&nbsp; std::cout &lt;&lt; &quot; Translation X = &quot; &lt;&lt; finalTranslationX&nbsp; &lt;&lt; std::endl;<br>&nbsp; std::cout &lt;&lt; &quot; Translation Y = &quot; &lt;&lt; finalTranslationY&nbsp; &lt;&lt; std::endl;<br>&nbsp; std::cout &lt;&lt; &quot; Iterations&nbsp;&nbsp;&nbsp; = &quot; &lt;&lt; numberOfIterations &lt;&lt; std::endl;
<br>&nbsp; std::cout &lt;&lt; &quot; Metric value&nbsp; = &quot; &lt;&lt; bestValue&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &lt;&lt; std::endl;<br>&nbsp; <br>&nbsp; vnl_matrix&lt;double&gt; p(2, 2);<br>&nbsp; p[0][0] = (double) finalParameters[0];<br>&nbsp; p[0][1] = (double) finalParameters[1];
<br>&nbsp; p[1][0] = (double) finalParameters[2];<br>&nbsp; p[1][1] = (double) finalParameters[3];<br>&nbsp; vnl_svd&lt;double&gt; svd(p);<br>&nbsp; vnl_matrix&lt;double&gt; r(2, 2);<br>&nbsp; r = svd.U() * vnl_transpose(svd.V());<br>&nbsp; double angle = asin(r[1][0]);
<br>&nbsp; <br>&nbsp; std::cout &lt;&lt; &quot; Scale 1&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = &quot; &lt;&lt; svd.W(0)&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &lt;&lt; std::endl;<br>&nbsp; std::cout &lt;&lt; &quot; Scale 2&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; = &quot; &lt;&lt; svd.W(1)&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &lt;&lt; std::endl;
<br>&nbsp; std::cout &lt;&lt; &quot; Angle (degrees) = &quot; &lt;&lt; angle * 45.0 / atan(1.0) &lt;&lt; std::endl;<br>&nbsp; <br><br>&nbsp; typedef itk::ResampleImageFilter&lt; <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; MovingImageType, <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; FixedImageType &gt;&nbsp;&nbsp;&nbsp; ResampleFilterType;
<br><br>&nbsp; TransformType::Pointer finalTransform = TransformType::New();<br><br>&nbsp; finalTransform-&gt;SetCenter( transform-&gt;GetCenter() );<br>&nbsp; finalTransform-&gt;SetParameters( finalParameters );<br><br>&nbsp; ResampleFilterType::Pointer resampler = ResampleFilterType::New();
<br><br>&nbsp; resampler-&gt;SetTransform( finalTransform );<br>&nbsp; resampler-&gt;SetInput( movingImageReader-&gt;GetOutput() );<br><br>&nbsp; FixedImageType::Pointer fixedImage = fixedImageReader-&gt;GetOutput();<br><br>&nbsp; resampler-&gt;SetSize(&nbsp;&nbsp;&nbsp; fixedImage-&gt;GetLargestPossibleRegion().GetSize() );
<br>&nbsp; resampler-&gt;SetOutputOrigin(&nbsp; fixedImage-&gt;GetOrigin() );<br>&nbsp; resampler-&gt;SetOutputSpacing( fixedImage-&gt;GetSpacing() );<br>&nbsp; resampler-&gt;SetDefaultPixelValue( 100 );<br>&nbsp; <br>&nbsp; typedef&nbsp; unsigned char&nbsp; OutputPixelType;
<br><br>&nbsp; typedef itk::Image&lt; OutputPixelType, Dimension &gt; OutputImageType;<br>&nbsp; <br>&nbsp; typedef itk::CastImageFilter&lt; <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; FixedImageType,<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OutputImageType &gt; CastFilterType;
<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <br>&nbsp; typedef itk::ImageFileWriter&lt; OutputImageType &gt;&nbsp; WriterType;<br><br><br>&nbsp; WriterType::Pointer&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; writer =&nbsp; WriterType::New();<br>&nbsp; CastFilterType::Pointer&nbsp; caster =&nbsp; CastFilterType::New();
<br><br><br>&nbsp; writer-&gt;SetFileName( argv[3] );<br><br><br>&nbsp; caster-&gt;SetInput( resampler-&gt;GetOutput() );<br>&nbsp; writer-&gt;SetInput( caster-&gt;GetOutput()&nbsp;&nbsp; );<br>&nbsp; writer-&gt;Update();<br><br><br>&nbsp; typedef itk::SubtractImageFilter&lt; 
<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; FixedImageType, <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; FixedImageType, <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; FixedImageType &gt; DifferenceFilterType;<br><br>&nbsp; DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
<br><br>&nbsp; difference-&gt;SetInput1( fixedImageReader-&gt;GetOutput() );<br>&nbsp; difference-&gt;SetInput2( resampler-&gt;GetOutput() );<br><br>&nbsp; WriterType::Pointer writer2 = WriterType::New();<br>&nbsp; <br>&nbsp; typedef itk::RescaleIntensityImageFilter&lt; 
<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; FixedImageType, <br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; OutputImageType &gt;&nbsp;&nbsp; RescalerType;<br><br>&nbsp; RescalerType::Pointer intensityRescaler = RescalerType::New();<br><br>&nbsp; intensityRescaler-&gt;SetInput( difference-&gt;GetOutput() );
<br>&nbsp; intensityRescaler-&gt;SetOutputMinimum(&nbsp;&nbsp; 0 );<br>&nbsp; intensityRescaler-&gt;SetOutputMaximum( 255 );<br>&nbsp; <br>&nbsp; writer2-&gt;SetInput( intensityRescaler-&gt;GetOutput() );&nbsp; <br>&nbsp; resampler-&gt;SetDefaultPixelValue( 1 );
<br>&nbsp;<br><br>&nbsp; if( argc &gt; 5 )<br>&nbsp;&nbsp;&nbsp; {<br>&nbsp;&nbsp;&nbsp; writer2-&gt;SetFileName( argv[5] );<br>&nbsp;&nbsp;&nbsp; writer2-&gt;Update();<br>&nbsp;&nbsp;&nbsp; }<br><br><br>&nbsp; typedef itk::IdentityTransform&lt; double, Dimension &gt; IdentityTransformType;<br>&nbsp; IdentityTransformType::Pointer identity = IdentityTransformType::New();
<br><br>&nbsp; if( argc &gt; 4 )<br>&nbsp;&nbsp;&nbsp; {<br>&nbsp;&nbsp;&nbsp; resampler-&gt;SetTransform( identity );<br>&nbsp;&nbsp;&nbsp; writer2-&gt;SetFileName( argv[4] );<br>&nbsp;&nbsp;&nbsp; writer2-&gt;Update();<br>&nbsp;&nbsp;&nbsp; }<br><br><br>&nbsp; return 0;<br>}<br><br>==================================================================
<br>