No subject

to register both images satisfactorily,... is that correct ?


    Please let us know,


         Thanks


             Luis



---------------------------------------------------------------------------=
------------------
On Sat, Nov 28, 2009 at 12:00 PM, Auron Ford Huton <auron at hotmail.es> wrote=
:
>
> Hi Luis and ITK!!
> First of all, really thank you for the time you expend answering me.
> About the things you told me:
>
> 1) You are right, If I don=B4t want the center to be calculated, is bette=
r to
> used Rigid2DTramsform, I have already changed it in my code, thank you.
>
> 2) What can I do to fixed the problem with the orientation? As I said, th=
e
> only way I can figure it out is by =A0"optimizerScales[0] =3D -1.0;" but =
like
> you said I sould never do that, is there other choices?
>
> 3) My images are about 40 consecutive slices from a microscopic tissue (t=
he
> ones I added as a screenshot were two of them) and I have to align them a=
ll
> so I can make after a 3D representation with VTK. So it doesn=B4t really
> matter wich one is the fixed or wich one is the moving. Why you asked me
> about it? Is there something that you saw that make you think I should us=
e
> them backwards? (I=B4m attaching a bigger screenshot, just in case it hel=
ps).
>
> 4) The rotation I need is about 50 degrees (not 90 as you told me) using =
as
> center of rotation the center of the image. By the way, the translation I
> calculated the images need is about -973 in the X axis and 2436 in the Y
> axis.
>
> Thanks for everything and sorry if my questions are a bit simple, but I=
=B4m a
> kind of a newby in ITK and Image Proccesing.
>
> P.D.: In the screenshot I attached, the left image is the fixed one and t=
he
> right image is the moving one.
>
> I also add below the code of my main method:
>
>
> int main( int argc, char *argv[] ){
> =A0=A0if( argc < 4 )=A0{
> =A0=A0 =A0std::cerr << "Missing Parameters " << std::endl;
> =A0=A0 =A0std::cerr << "Usage: " << argv[0];
> =A0=A0 =A0std::cerr << " fixedImageFile =A0movingImageFile ";
> =A0=A0 =A0std::cerr << " outputImagefile =A0[differenceOutputfile] ";
> =A0=A0 =A0std::cerr << " [differenceBeforeRegistration] "<< std::endl;
> =A0=A0 =A0return EXIT_FAILURE;
> =A0=A0 =A0}
>
>
> =A0=A0const =A0 =A0unsigned int =A0 =A0Dimension =3D 2;
> =A0=A0typedef =A0unsigned char =A0 =A0 =A0 =A0 =A0 PixelType;
> =A0=A0typedef itk::OrientedImage< PixelType, Dimension > =A0FixedImageTyp=
e;
> =A0=A0typedef itk::OrientedImage< PixelType, Dimension > =A0MovingImageTy=
pe;
> =A0=A0typedef itk::Rigid2DTransform< double > TransformType;
>
>
> =A0=A0typedef itk::RegularStepGradientDescentOptimizer =A0 =A0 =A0 Optimi=
zerType;
> =A0=A0typedef itk::MeanSquaresImageToImageMetric<
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0FixedImageType,
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0MovingImageType > =A0 =A0MetricType;
> =A0=A0typedef itk:: LinearInterpolateImageFunction<
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0MovingImageType,
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0double =A0 =A0 =A0 =A0 =A0> =A0 =A0InterpolatorType;
> =A0=A0typedef itk::ImageRegistrationMethod<
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0FixedImageType,
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0MovingImageType > =A0 =A0RegistrationType;
> =A0=A0MetricType::Pointer =A0 =A0 =A0 =A0 metric =A0 =A0 =A0 =A0=3D Metri=
cType::New();
> =A0=A0OptimizerType::Pointer =A0 =A0 =A0optimizer =A0 =A0 =3D OptimizerTy=
pe::New();
> =A0=A0InterpolatorType::Pointer =A0 interpolator =A0=3D InterpolatorType:=
:New();
> =A0=A0RegistrationType::Pointer =A0 registration =A0=3D RegistrationType:=
:New();
>
> =A0=A0registration->SetMetric( =A0 =A0 =A0 =A0metric =A0 =A0 =A0 =A0);
> =A0=A0registration->SetOptimizer( =A0 =A0 optimizer =A0 =A0 );
> =A0=A0registration->SetInterpolator( =A0interpolator =A0);
>
> =A0=A0TransformType::Pointer =A0transform =3D TransformType::New();
> =A0=A0registration->SetTransform( transform );
>
> =A0=A0typedef itk::ImageFileReader< FixedImageType =A0> FixedImageReaderT=
ype;
> =A0=A0typedef itk::ImageFileReader< MovingImageType > MovingImageReaderTy=
pe;
> =A0=A0FixedImageReaderType::Pointer =A0fixedImageReader =A0=3D
> FixedImageReaderType::New();
> =A0=A0MovingImageReaderType::Pointer movingImageReader =3D
> MovingImageReaderType::New();
> =A0=A0fixedImageReader->SetFileName( =A0argv[1] );
> =A0=A0movingImageReader->SetFileName( argv[2] );
>
> =A0=A0registration->SetFixedImage( =A0 =A0fixedImageReader->GetOutput() =
=A0 =A0);
> =A0=A0registration->SetMovingImage( =A0 movingImageReader->GetOutput() =
=A0 );
> =A0=A0fixedImageReader->Update();
> =A0=A0registration->SetFixedImageRegion(
> =A0=A0 =A0 fixedImageReader->GetOutput()->GetBufferedRegion() );
> =A0=A0typedef itk::CenteredTransformInitializer<
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0TransformType,
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0FixedImageType,
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =
=A0MovingImageType >
> =A0TransformInitializerType;
> =A0=A0TransformInitializerType::Pointer initializer =3D
> TransformInitializerType::New();
>
> =A0=A0initializer->SetTransform( =A0 transform );
> =A0=A0initializer->SetFixedImage( =A0fixedImageReader->GetOutput() );
> =A0=A0initializer->SetMovingImage( movingImageReader->GetOutput() );
>
> =A0=A0initializer->MomentsOn();
> =A0=A0//initializer->GeometryOn();
>
> =A0=A0initializer->InitializeTransform();
>
> =A0=A0transform->SetAngle( 0.0 );//0.0
>
> =A0=A0const FixedImageType::SpacingType&
> =A0=A0 =A0spacing =3D fixedImageReader->GetOutput()->GetSpacing();
> =A0=A0const FixedImageType::PointType&
> =A0=A0 =A0origin =A0=3D fixedImageReader->GetOutput()->GetOrigin();
> =A0=A0FixedImageType::SizeType size =3D
> =A0=A0 =A0 =A0fixedImageReader->GetOutput()->GetLargestPossibleRegion().G=
etSize();
> =A0=A0TransformType::InputPointType rotationCenter;
> =A0=A0rotationCenter[0] =3D origin[0] + spacing[0] * size[0] / 2.0;
> =A0=A0rotationCenter[1] =3D origin[1] + spacing[1] * size[1] / 2.0;
> =A0=A0std::cout<< "RotationCeterX" << origin[0] + spacing[0] * size[0] / =
2.0
> <<std::endl;
> =A0=A0std::cout<< "RotationCeterY" << origin[1] + spacing[1] * size[1] / =
2.0
> <<std::endl;
> =A0=A0transform->SetCenter( rotationCenter );
>
>
> =A0=A0registration->SetInitialTransformParameters( transform->GetParamete=
rs() );
>
> =A0=A0typedef OptimizerType::ScalesType =A0 =A0 =A0 OptimizerScalesType;
> =A0=A0OptimizerScalesType optimizerScales( transform->GetNumberOfParamete=
rs() );
> =A0=A0const double translationScale =3D 1.0 / 5000.0;
> =A0=A0std::cout<< "translationScale" << translationScale <<std::endl;
> =A0=A0optimizerScales[0] =3D 1.0;
> =A0=A0optimizerScales[1] =3D translationScale;
> =A0=A0optimizerScales[2] =3D translationScale;
>
> =A0=A0optimizer->SetScales( optimizerScales );
> =A0=A0optimizer->SetMaximumStepLength(0.1);
> =A0=A0optimizer->SetMinimumStepLength(0.001 );
> =A0=A0optimizer->SetNumberOfIterations(200);
> =A0=A0//optimizer->SetRelaxationFactor( 0.9 );
>
> =A0=A0CommandIterationUpdate::Pointer observer =3D CommandIterationUpdate=
::New();
> =A0=A0optimizer->AddObserver( itk::IterationEvent(), observer );
> =A0=A0try
> =A0=A0 =A0{
> =A0=A0 =A0registration->StartRegistration();
> =A0=A0 =A0}
> =A0=A0catch( itk::ExceptionObject & err )
> =A0=A0 =A0{
> =A0=A0 =A0std::cerr << "ExceptionObject caught !" << std::endl;
> =A0=A0 =A0std::cerr << err << std::endl;
> =A0=A0 =A0return EXIT_FAILURE;
> =A0=A0 =A0}
>
> =A0=A0OptimizerType::ParametersType finalParameters =3D
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0registration->GetLastTransformP=
arameters();
>
>
> =A0=A0const double finalAngle =A0 =A0 =A0 =A0 =A0 =3D finalParameters[0];
> =A0=A0const double finalTranslationX =A0 =A0=3D finalParameters[1];
> =A0=A0const double finalTranslationY =A0 =A0=3D finalParameters[2];
> =A0=A0const unsigned int numberOfIterations =3D optimizer->GetCurrentIter=
ation();
> =A0=A0const double bestValue =3D optimizer->GetValue();
> =A0=A0// Print out results
> =A0=A0//
> =A0=A0const double finalAngleInDegrees =3D finalAngle * 45.0 / atan(1.0);
> =A0=A0std::cout << "Result =3D " << std::endl;
> =A0=A0std::cout << " Angle (radians) " << finalAngle =A0<< std::endl;
> =A0=A0std::cout << " Angle (degrees) " << finalAngleInDegrees =A0<< std::=
endl;
> =A0=A0std::cout << " Center X =A0 =A0 =A0=3D " << rotationCenter[0] =A0<<=
 std::endl;
> =A0=A0std::cout << " Center Y =A0 =A0 =A0=3D " << rotationCenter[1] =A0<<=
 std::endl;
> =A0=A0std::cout << " Translation X =3D " << finalTranslationX =A0<< std::=
endl;
> =A0=A0std::cout << " Translation Y =3D " << finalTranslationY =A0<< std::=
endl;
> =A0=A0std::cout << " Iterations =A0 =A0=3D " << numberOfIterations << std=
::endl;
> =A0=A0std::cout << " Metric value =A0=3D " << bestValue =A0 =A0 =A0 =A0 =
=A0<< std::endl;
> =A0=A0std::cout << " Stop Condition =A0=3D " << optimizer->GetStopConditi=
on() <<
> std::endl;
>
> =A0=A0transform->SetParameters( finalParameters );
> =A0=A0TransformType::MatrixType matrix =3D transform->GetRotationMatrix()=
;
> =A0=A0TransformType::OffsetType offset =3D transform->GetOffset();
> =A0=A0std::cout << "Matrix =3D " << std::endl << matrix << std::endl;
> =A0=A0std::cout << "Offset =3D " << std::endl << offset << std::endl;
>
>
> =A0=A0typedef itk::ResampleImageFilter<
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0MovingImageType=
,
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0FixedImageType =
> =A0 =A0ResampleFilterType;
> =A0=A0TransformType::Pointer finalTransform =3D TransformType::New();
> =A0=A0finalTransform->SetParameters( finalParameters );
> =A0=A0ResampleFilterType::Pointer resample =3D ResampleFilterType::New();
> =A0=A0resample->SetTransform( finalTransform );
> =A0=A0resample->SetInput( movingImageReader->GetOutput() );
> =A0=A0FixedImageType::Pointer fixedImage =3D fixedImageReader->GetOutput(=
);
> =A0=A0resample->SetSize( =A0 =A0fixedImage->GetLargestPossibleRegion().Ge=
tSize() );
> =A0=A0resample->SetOutputOrigin( =A0fixedImage->GetOrigin() );
> =A0=A0resample->SetOutputSpacing( fixedImage->GetSpacing() );
> =A0=A0resample->SetOutputDirection( fixedImage->GetDirection() );
> =A0=A0resample->SetDefaultPixelValue( 100 );
>
> =A0=A0typedef =A0unsigned char =A0OutputPixelType;
> =A0=A0typedef itk::OrientedImage< OutputPixelType, Dimension > OutputImag=
eType;
>
> =A0=A0typedef itk::CastImageFilter<
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0FixedImageType,
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0OutputImageType > CastF=
ilterType;
>
> =A0=A0typedef itk::ImageFileWriter< OutputImageType > =A0WriterType;
>
> =A0=A0WriterType::Pointer =A0 =A0 =A0writer =3D =A0WriterType::New();
> =A0=A0CastFilterType::Pointer =A0caster =3D =A0CastFilterType::New();
>
> =A0=A0writer->SetFileName( argv[3] );
>
> =A0=A0caster->SetInput( resample->GetOutput() );
> =A0=A0writer->SetInput( caster->GetOutput() =A0 );
> =A0=A0writer->Update();
>
> =A0=A0typedef itk::SquaredDifferenceImageFilter<
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0Fix=
edImageType,
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0Fix=
edImageType,
> =A0=A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0Out=
putImageType > DifferenceFilterType;
> =A0=A0DifferenceFilterType::Pointer difference =3D DifferenceFilterType::=
New();
> =A0=A0WriterType::Pointer writer2 =3D WriterType::New();
> =A0=A0writer2->SetInput( difference->GetOutput() );
>
> =A0=A0if( argc >=3D 5 )
> =A0=A0 =A0{
> =A0=A0 =A0difference->SetInput1( fixedImageReader->GetOutput() );
> =A0=A0 =A0difference->SetInput2( resample->GetOutput() );
> =A0=A0 =A0writer2->SetFileName( argv[4] );
> =A0=A0 =A0writer2->Update();
> =A0=A0 =A0}
>
> =A0=A0if( argc >=3D 6 )
> =A0=A0 =A0{
> =A0=A0 =A0writer2->SetFileName( argv[5] );
> =A0=A0 =A0difference->SetInput1( fixedImageReader->GetOutput() );
> =A0=A0 =A0difference->SetInput2( movingImageReader->GetOutput() );
> =A0=A0 =A0writer2->Update();
> =A0=A0 =A0}
> =A0=A0return EXIT_SUCCESS;
> }
>
>
>
> Thank you in advance.
>
> Auron.
>> Date: Thu, 26 Nov 2009 14:34:58 -0500
>> Subject: Re: [Insight-users] Confused abour Optimizer Scales
>> From: luis.ibanez at kitware.com
>> To: auron at hotmail.es
>> CC: insight-users at itk.org
>>
>> Hi Auron,
>>
>>
>> The class
>>
>> itkCenteredRigid2DTransform
>>
>> has five parameters, that are:
>>
>> 1) Angle
>> 2) Center of Rotation (X coordinate)
>> 3) Center of Rotation (X coordinate)
>> 4) Translation (X coordinate)
>> 5) Translation (Y coordinate)
>>
>> Several things are wrong with the current optimizer scales:
>>
>> > optimizerScales[0] =3D -1.0; // If I don=B4t use a negative value, the
>> > results
>> > turn negative.
>> > optimizerScales[1] =3D 1.0; // My center of rotation is the center of =
the
>> > image and I don=B4t want it to vary.
>> > optimizerScales[2] =3D 1.0;
>> > optimizerScales[3] =3D translationScale;
>> > optimizerScales[4] =3D translationScale;
>>
>>
>> A) the first one should NOT be negative.
>> Something else is going on with the orientation of your image.
>> but... in any case, the optimizer scales should NEVER be negative.
>>
>> B) The second and third values correspond to the center of rotation,
>> but you don't want to change this center. My suggestion will be
>> to put in here, very small numbers. For example 1e-6.
>>
>>
>> Also, since you don't want to modify the rotation center,
>> you may want to simply not use this Transform, and instead
>> use the class
>> itkRigid2DTransform.h
>>
>> that also offers a Center of rotation, but does not include it
>> as part of the parameters to be optimized.
>>
>>
>> Finally,
>> Thanks for adding the screenshots of your images.
>>
>> Are they really the Fixed image and Moving image ?
>>
>> If so,
>> please note that there NO WAY that the registration process
>> could compensate for their current rotation of 90 degrees :-)
>>
>> You have to give it a hand at the registration process and initialize
>> that angle with PI / 2. Otherwise, The Optimizer simply will not walk
>> that far in the space of the Transform parameters. Usually a registratio=
n
>> process will be able to deal with rotations below 20 degrees. (unless,
>> of course, you really customize them for a large rotation setting).
>>
>>
>>
>> Regards,
>>
>>
>> Luis
>>
>>
>>
>
>
> ________________________________
> =A1Nuevo Canal Mujer! Moda, belleza, sexo, dietas, embarazo. m=E1s f=E1ci=
l y a tu
> alcance. Si quieres estar a la =FAltima, no puedes perd=E9rtelo.