<div dir="ltr"><div><div><div>Hi Jérôme,<br><br></div>Have you tried changing the optimizer scales? <br>That parameter changes the sensitivity of the optimization process between the translation and rotation of the transform, it seams to me that the value 1.0/1e7 is probably too small and the optimizer is not sensitive enough to rotations.<br>
<br></div>For test purposes you can also set the mutual information metric to compute histograms on the whole image, that may give you a much detailed metric map for the optimizer as well.<br><br></div><div><span style="color:rgb(0,0,0)">metric</span><span style="color:rgb(0,0,0)">-></span>SetNumberOfSpatialSamples<span style="color:rgb(0,0,0)">(</span><span style="color:rgb(192,192,192)"> <some value></span><span style="color:rgb(0,0,128)"></span><span style="color:rgb(192,192,192)"> </span><span style="color:rgb(0,0,0)">);</span></div>
<div><br></div>Hope that helps,<br><div class="gmail_extra"><br clear="all"><div><div dir="ltr">Nicolás Gallego-Ortiz<br>Université catholique de Louvain, Belgium<br></div></div>
<br><br><div class="gmail_quote">2014-07-03 1:29 GMT+02:00 Jerome Plumat <span dir="ltr"><<a href="mailto:j.plumat@auckland.ac.nz" target="_blank">j.plumat@auckland.ac.nz</a>></span>:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Hi everyone,<br>
I'm facing to a problem with 3D affine multi resolution registration.<br>
The pipeline registers two 3D MRI head and necks volumes (600x400x40). These are scans of the same patient at two times. Due to the protocol, volumes are similar expect for MRI artifacts. In fine the algorithm would be used to registered an atlas defined on the fixed image on moving image.<br>
<br>
The main elements of the pipeline are listed below. It's implemented with itk4 and is mostly based on MutliResImageRegistration2.<u></u>cxx, I use the Mutual Information metric, linear interpolator, the optimizer is RegularStepGradientDescentOpti<u></u>mizer and CenteredTransformInitializer.<br>
<br>
The registration is performed without any error and the translation is correctly found but the algorithm can't provide any rotation and I can't figure why.<br>
The final results are<br>
versor X = [1.00007, 0.000485777, 1.5912e-05]<br>
versor Y = [-0.000126483, 0.999922, 1.34996e-05]<br>
versor Z = [2.35376e-05, -8.28481e-05, 0.999986]<br>
Translation X = 1.35204<br>
Translation Y = 0.647584<br>
Translation Z = -0.132708<br>
Iterations = 96<br>
Metric value = -0.738412<br>
<br>
Following <a href="http://www.cmake.org/pipermail/insight-users/2006-August/019025.html" target="_blank">http://www.cmake.org/<u></u>pipermail/insight-users/2006-<u></u>August/019025.html</a> the final rotation is 0.0003 along axis [-0.155456, -0.0122989, -0.987766].<br>
Moreover, visual validation clearly highlight the fact the rotation is not sufficient.<br>
<br>
Hope any one has advices.<br>
Thanks in advances.<br>
<br>
Jerome.<br>
<br>
<br>
Part of code:<br>
<br>
....<br>
<br>
// transform<br>
typedef itk::AffineTransform< double, Dimension> TransformType;<br>
typedef itk::<u></u>CenteredTransformInitializer< TransformType, FixedImageType, MovingImageType > TransformInitializerType;<br>
// optimizer<br>
typedef itk::<u></u>RegularStepGradientDescentOpti<u></u>mizer OptimizerType;<br>
<br>
// metric<br>
typedef itk::<u></u>MattesMutualInformationImageTo<u></u>ImageMetric<<br>
MovingImageType,<br>
MovingImageType > MetricType;<br>
// interpolator<br>
typedef itk:: LinearInterpolateImageFunction<u></u><<br>
MovingImageType,<br>
double > InterpolatorType;<br>
<br>
// registration<br>
typedef itk::<u></u>MultiResolutionImageRegistrati<u></u>onMethod<<br>
FixedImageType,<br>
MovingImageType > RegistrationType;<br>
<br>
<br>
typedef itk::<u></u>MultiResolutionPyramidImageFil<u></u>ter<<br>
FixedImageType,<br>
FixedImageType > FixedImagePyramidType;<br>
typedef itk::<u></u>MultiResolutionPyramidImageFil<u></u>ter<<br>
MovingImageType,<br>
MovingImageType > MovingImagePyramidType;<br>
<br>
<br>
...<br>
<br>
metric-><u></u>SetNumberOfHistogramBins( 256 );<br>
metric-><u></u>SetNumberOfSpatialSamples( 50000 );<br>
metric->ReinitializeSeed<a href="tel:%28%2076926294" value="+5776926294" target="_blank">( 76926294</a> );<br>
<br>
....<br>
<br>
// connect the inputs<br>
registration->SetFixedImage( fixedImage );<br>
registration->SetMovingImage( movingImage );<br>
<br>
<br>
// set the regions<br>
registration-><u></u>SetFixedImageRegion( fixedImage-><u></u>GetLargestPossibleRegion() );<br>
<br>
// configure the initializer<br>
initializer->SetTransform( transform );<br>
initializer->SetFixedImage( fixedImage );<br>
initializer->SetMovingImage( movingImage );<br>
initializer->MomentsOn();<br>
initializer-><u></u>InitializeTransform();<br>
<br>
registration-><u></u>SetInitialTransformParameters( transform->GetParameters() );<br>
<br>
// optimizer scales<br>
typedef OptimizerType::ScalesType OptimizerScalesType;<br>
OptimizerScalesType optimizerScales( transform-><u></u>GetNumberOfParameters() );<br>
const double momentScale = 1;<br>
const double translationScale = 1.0 / 1e7;<br>
<br>
//number of dimensions: N*(N+1)<br>
for(unsigned int i=0; i<(Dimension*Dimension); i++ )<br>
{<br>
optimizerScales[i] = momentScale;<br>
}<br>
<br>
for( unsigned int i=(Dimension*Dimension);<br>
i<((Dimension+1)*Dimension); i++ )<br>
{<br>
optimizerScales[i] = translationScale;<br>
}<br>
optimizer->SetScales( optimizerScales );<br>
<br>
<br>
// performed registration<br>
...<br>
<br>
//<br>
// deformed the moving image and save it<br>
//<br>
typedef itk::ResampleImageFilter<<br>
MovingImageType,<br>
FixedImageType > ResampleFilterType;<br>
typedef itk::CastImageFilter<<br>
MovingImageType,<br>
OutputImageType > CastFilterType;<br>
typedef itk::ImageFileWriter< OutputImageType > WriterType;<br>
<br>
ResampleFilterType::Pointer resampler = ResampleFilterType::New();<br>
WriterType::Pointer writer = WriterType::New();<br>
CastFilterType::Pointer caster = CastFilterType::New();<br>
<br>
TransformType::Pointer finalTransform = TransformType::New();<br>
<br>
finalTransform->SetParameters(<br>
registration-><u></u>GetLastTransformParameters() );<br>
finalTransform-><u></u>SetFixedParameters(<br>
transform->GetFixedParameters(<u></u>) );<br>
<br>
// set transform and initial moving image<br>
resampler->SetTransform( finalTransform );<br>
resampler->SetInput( movingImage );<br>
<br>
resampler->SetSize( fixedImage-><u></u>GetLargestPossibleRegion().<u></u>GetSize() );<br>
resampler->SetOutputOrigin( fixedImage->GetOrigin() );<br>
resampler->SetOutputSpacing( fixedImage->GetSpacing() );<br>
resampler->SetOutputDirection( fixedImage->GetDirection() );<br>
resampler-><u></u>SetDefaultPixelValue( defaultOutputPixel );<br>
<br>
writer->SetFileName( argv[argvOutputImageFile] );<br>
caster->SetInput( resampler->GetOutput() );<br>
writer->SetInput( caster->GetOutput() );<br>
<br>
...<br>
<br>
-- <br>
<br>
Jerome<br>
-------<br>
School of Medical Sciences<br>
University of Auckland<br>
<br>
If I am not for myself, who will be for me? And if I am only for myself, then what an I? And if not now when? – Hillel HaZaken<br>
<br>
______________________________<u></u>_______<br>
Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
<br>
Visit other Kitware open-source projects at<br>
<a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/<u></u>opensource/opensource.html</a><br>
<br>
Kitware offers ITK Training Courses, for more information visit:<br>
<a href="http://www.kitware.com/products/protraining.php" target="_blank">http://www.kitware.com/<u></u>products/protraining.php</a><br>
<br>
Please keep messages on-topic and check the ITK FAQ at:<br>
<a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_<u></u>FAQ</a><br>
<br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://public.kitware.com/mailman/listinfo/insight-users" target="_blank">http://public.kitware.com/<u></u>mailman/listinfo/insight-users</a><br>
______________________________<u></u>_________________<br>
Community mailing list<br>
<a href="mailto:Community@itk.org" target="_blank">Community@itk.org</a><br>
<a href="http://public.kitware.com/mailman/listinfo/community" target="_blank">http://public.kitware.com/<u></u>mailman/listinfo/community</a><br>
</blockquote></div><br></div></div>