[Insight-users] Registration of US and CT images
Anders Almås
andersal@stud.ntnu.no
Wed May 12 10:05:50 EDT 2004
Hi Luis!
Thanks for the help about US - CT registration. What you wrote was very
helpful!!
And about the step length in the optimizer which I forgot to send along
the last mail, those are just what you where writing. I am using a
pyramide with 3-4 levels and at the starting level the max step length
is 1 and min step length is 0.1 . At each new level the max step is
decreased with the resulting step length in the last level and the new
min level is 1/10 of the last levels min length.
Regards
Anders
> Subject: Re: [Insight-users] Registration of US and CT images
>
>
> Hi Anders
>
> The easy way to initialize the VersorRigid3DTransform
> is to first call SetIdentity() on it, and then use
> the CenteredTransformInitializer in order to set the
> initial translation.
>
> The Versor is not limited to manage rotations in one
> direction. It will represent a rotation around an
> arbitrary axis in 3D space.
>
> Note that the VersorRigid3D transform can only represent
> Rigid transforms not generic Affine transforms. If you
> applied an Affine transform to your dataset and this
> Affine transform included anisotropic scaling or shearing,
> the VersorRigid3DTransform *will not* be able to correct
> for those deformations. You are free to use an Affine
> transform for initializing the dataset as long as you
> only set rotations and translation in that transform.
>
> Actually, if you want to configure a test, the easy way
> to go is to miss-register the image on purpose using
> a VersorRigid3DTransform. Then you perform the registration
> and verify if the solution is the inverse of the transform
> that you manually applied at the beginning.
>
>
> ---
>
>
> The optimzerScales are the usual suspect in any
> problem related to Transforms combining translations,
> rotations and scaling. Your values seems to be ok for
> the scale of the image though.
>
> The next usual suspect is the step length of the optimizer.
> I couldn't find in your code the setup for this parameters.
>
> Did I missed it ?
> or is that you omitted to set this parameter ?
>
> The step length is critical in the optimization process.
> You should start with very small values, for example
> 0.01, and only increase them if you see that the transform
> parameters move too slowly.
>
> The best way for you to get a grip on the appropriate
> parameters for the optimizer is to manually subsample the
> images to the lowest resolution level of the pyramid, and
> feed those images in a single-resolution level program.
> Fine tune the parameters of the optimizer for that level
> until you get convergence.
>
> Then you can move to the next resolution level, and so on...
>
> Note that *there is no reason* to expect that the optimization
> parameters that were useful at one resolution level will
> also be appropriate for another resolution level. They are
> a good starting point though.
>
> Once you figure out the parameters for every resolution level
> you can come back to the multi-resolution level code and
> make the Command/Observer set up the appropriate parameters
> for each level as part of its response to the iteration (the
> change of resolution level).
>
>
>
> Regards,
>
>
>
> Luis
>
>
>
> -----------------
> Anders Almås wrote:
>
> > Hi!
> >
> > I am performing registration of US images and CT images in 3D. So
> far I
> > am only using Mattes mutual information metric. Now my ragistration
> is
> > only performed on a perfect aligned set of US and CT images. What I
> do
> > is to initially rotate and translate the CT image and try to
> registrate
> > it such that the translated/rotated image is aligned again. First
> i
> > used the Centered Affine transform, but that gave only results
> which
> > was sometimes ok for only rotation. Then when I only tested for
> > translations the translation transform worked quite good. At
> present
> > time I am testing the VersorRigid3DTransform. This one works ok at
> > translations and it seems that it manages rotations in one
> direction.
> > My question for VersorRigid3DTransform is how it should be
> initialized?
> > I send the code so that any one can comment it... Is it the
> > optimizerscales which are the problem. Or is it really that I
> perform
> > an affine transformation/rotation that makes it difficult to
> registrate
> > it back? I hope anyone can help me about this...
> > The size of the images is 256*256*176 with spacing 1.453,1.453 and
> 2
> >
> > Here follows the code for the registration:
> >
> > std::cout<<"Here starts the registration process..."<<std::endl;
> >
> >
> > //TranslationTransform
> > // typedef itk::TranslationTransform<double, ImageDimension>
> > TranslationType;
> > //Rigid3DTransform
> > typedef itk::VersorRigid3DTransform<double> TranslationType;
> > // typedef itk::CenteredAffineTransform< double, ImageDimension >
>
> > TranslationType;
> > //Transform filters
> > //Optimizer filters
> > typedef itk::VersorRigid3DTransformOptimizer OptimizerType;
> >
> > // typedef itk::RegularStepGradientDescentOptimizer
> > OptimizerType;
> >
> > //Image Function filters
> > typedef itk::LinearInterpolateImageFunction<ImageType, double >
> > InterpolatorType;
> >
> >
> > //Metric filter
> > typedef itk::MattesMutualInformationImageToImageMetric<
> ImageType,
> > ImageType > MetricType;
> > //typedef
> itk::MutualInformationImageToImageMetric<ImageType,ImageType
> >
> >>MetricType;
> >
> >
> > //RegistrationMethod
> > //typedef itk::ImageRegistrationMethod< ImageType, ImageType
> >
> > RegistrationType;
> > typedef itk::MultiResolutionImageRegistrationMethod<
> ImageType,
> > ImageType > RegistrationType;
> >
> > typedef
> itk::MultiResolutionPyramidImageFilter<ImageType,ImageType >
> > FixedImagePyramidType;
> > typedef
> itk::MultiResolutionPyramidImageFilter<ImageType,ImageType >
> > MovingImagePyramidType;
> >
> > typedef OptimizerType::ScalesType OptimizerScalesType;
> >
> > TranslationType::Pointer transform2 = TranslationType::New();
> > OptimizerType::Pointer optimizer =
> OptimizerType::New();
> > InterpolatorType::Pointer interpolator2 =
> InterpolatorType::New();
> > RegistrationType::Pointer registration =
> RegistrationType::New();
> > MetricType::Pointer metric = MetricType::New();
> >
> >
> > FixedImagePyramidType::Pointer fixedImagePyramid =
> > FixedImagePyramidType::New();
> > MovingImagePyramidType::Pointer movingImagePyramid =
> > MovingImagePyramidType::New();
> >
> > typedef itk::CenteredTransformInitializer< TranslationType,
> ImageType,
> > ImageType > TransformInitializerType;
> > TransformInitializerType::Pointer initializer =
> > TransformInitializerType::New();
> > initializer->SetTransform( transform2 );
> > initializer->SetFixedImage( range->GetOutput());
> > // initializer->SetFixedImage( gaussian->GetOutput());
> > initializer->SetMovingImage( filter->GetOutput() );
> > //initializer->SetMovingImage( ctimage );
> > initializer->MomentsOn();
> > initializer->InitializeTransform();
> >
> > typedef TranslationType::VersorType VersorType;
> > typedef VersorType::VectorType VectorType;
> >
> > VersorType rotation;
> > VectorType axis;
> >
> > axis[0] = 0.0;
> > axis[1] = 0.0;
> > axis[2] = 1.0;
> >
> > const double angle = 0;
> >
> > rotation.Set( axis, angle );
> >
> > transform2->SetRotation( rotation );
> > axis[0] = 1.0;
> > axis[1] = 0.0;
> > axis[2] = 0.0;
> >
> > rotation.Set( axis, angle );
> >
> > transform2->SetRotation( rotation );
> >
> > axis[0] = 0.0;
> > axis[1] = 1.0;
> > axis[2] = 0.0;
> >
> > // const double angle = 0;
> >
> > rotation.Set( axis, angle );
> >
> > transform2->SetRotation( rotation );
> > //transform2->SetIdentity();
> > registration->SetTransform(transform2);
> > registration->SetOptimizer( optimizer );
> > registration->SetInterpolator( interpolator2 );
> > registration->SetMetric( metric );
> > registration->SetFixedImagePyramid( fixedImagePyramid );
> > registration->SetMovingImagePyramid( movingImagePyramid );
> >
> > OptimizerScalesType optimizerScales(
> > transform2->GetNumberOfParameters() );
> >
> > std::cout<<"Number of parametres:
> > "<<transform2->GetNumberOfParameters()<<std::endl;
> >
> >
> > optimizerScales[0] = 1.0;
> > optimizerScales[1] = 1.0;
> > optimizerScales[2] = 1.0;
> > optimizerScales[3] = 1.0/1000;
> >
> > optimizerScales[4] = 1.0/1000;
> > optimizerScales[5] = 1.0/1000 ;
> >
> > optimizer->SetScales( optimizerScales );
> > optimizer->MaximizeOn();
> >
> > registration->SetFixedImage( range->GetOutput()); //US IMAGE
> EDGE
> > MAP
> > registration->SetMovingImage( filter->GetOutput() ); //
> CTimage
> > gaussian blurred translated and rotated...
> > registration->SetFixedImageRegion(
> > range->GetOutput()->GetBufferedRegion() );
> >
> > typedef RegistrationType::ParametersType ParametersType;
> >
> > /*
> > ParametersType initialParameters(
> transform->GetNumberOfParameters()
> > );
> >
> > initialParameters[0] = 0.0; // Initial offset in mm along X
> > initialParameters[1] = 0.0; // Initial offset in mm along Y
> > initialParameters[2] = 0.0; // Initial offset in mm along Z
> > */
> > registration->SetInitialTransformParameters(
> > transform2->GetParameters() );
> >
> > metric->SetNumberOfHistogramBins( par[7] ); // 20
> > metric->SetNumberOfSpatialSamples( par[8] ); //10000
> >
> > optimizer->SetNumberOfIterations( par[6] ); // Tested from 300 -
> 800
> > iterations
> >
> > CommandIterationUpdate::Pointer observer =
> > CommandIterationUpdate::New();
> > optimizer->AddObserver( itk::IterationEvent(), observer );
> >
> >
> >
> > typedef RegistrationInterfaceCommand<RegistrationType>
> CommandType;
> > CommandType::Pointer command = CommandType::New();
> > registration->AddObserver( itk::IterationEvent(), command );
> >
> > registration->SetNumberOfLevels( par[9] ); //Tested from 3-5
> levels..
> >
> > try
> > {
> > registration->StartRegistration();
> > }
> > catch( itk::ExceptionObject & err )
> > {
> > std::cout << "ExceptionObject caught !" << std::endl;
> > std::cout << err << std::endl;
> >
> > return input;
> > }
> >
> > ParametersType finalParameters =
> > registration->GetLastTransformParameters();
> >
> > unsigned int numberOfIterations =
> optimizer->GetCurrentIteration();
> >
> >
> > double bestValue = optimizer->GetValue();
> >
> >
> > // Print out results
> > //
> > std::cout << "Result = " << std::endl;
> > std::cout << " Iterations = " << numberOfIterations <<
> std::endl;
> > std::cout << " Metric value = " << bestValue <<
> std::endl;
> > //std::cout << " Stop Condition = " <<
> optimizer->GetStopCondition()
> > << std::endl;
> >
> > I really hope somebody can help me with some advice. Maybe is it
> that my
> > parameters is set wrong??
> >
> > Regards
> >
> > Anders
> >
> > _______________________________________________
> > Insight-users mailing list
> > Insight-users@itk.org
> > http://www.itk.org/mailman/listinfo/insight-users
> >
>
>
>
>
> --__--__--
>
> Message: 5
> Date: Tue, 11 May 2004 18:09:19 -0400
> From: Luis Ibanez <luis.ibanez@kitware.com>
> To: Vaijanath Narasinga Rao <vaiju@cse.iitb.ac.in>
> Cc: insight-users@itk.org
> Subject: Re: [Insight-users] Parallel Processing
>
>
> Hi Vaijanath,
>
> We will be happy to get your algorithms in ITK.
>
> Here is the typical process to follow in order
> to include code in ITK:
>
>
> 1) Check that your filters are using the ITK
> coding style. you will find a PDF document
> describing the style under
>
> Insight/Documentation/Style.pdf
>
>
> 2) Make sure that your code compiles and runs.
>
> It may sound obvious,... but there are software
> developers out there who tend to forget
> these little details :-)
>
> Provide a reasonable test for every class
> in your filters. Test shouldn't be too
> intensive since we run them in a continuous
> basis. Typically something in the range of
> seconds would be ok.
>
>
> 3) About how to invoke a mpicc compiler, you
> can simply configure your CMakeLists.txt
> file in order to use a different compiler.
> The CMake CUSTOM_COMMAND is a flexible way
> to go about it.
>
> You might want to look at the files
>
> CMake/Modules/FindMPI.cmake
>
> and to the CMakeLists.txt file in VTK under
>
> VTK/Parallel
>
> in order to get a felling on how to configure
> a project that uses MPI. VTK (as opposed to ITK)
> has support for MPI, therefore you will find
> a lot of useful information regarding the required
> CMake configuration.
>
>
>
> 4) The full documentation of ITK filters can be
> found at
>
> http://www.itk.org/Insight/Doxygen/html/classes.html
>
> Note that parallelism in ITK is implemented only
> for shared memory, not for distributed systems.
>
> The filters that support parallel execution (in the
> shared memory setup) are listed under:
>
>
>
http://www.itk.org/Insight/Doxygen/html/group__MultithreadingGroup.html
>
>
>
>
>
>
> Regards,
>
>
>
> Luis
>
>
> -----------------------------
> Vaijanath Narasinga Rao wrote:
>
> > Hi To all,
> >
> > I wanted to know few things, I haev written some parallel
> algorithms for
> > Edge detection, Watermarking Segmentation, Thresholding techniques.
> I
> > wanted to know that how can i add them to ITK library. Also how can
> i
> > chaneg make files to compile them via mpicc and not the normal C
> compiler.
> >
> > Also one more thing, Where can I get the list of all filters that
> > (frequency as well as domain along with there documentation) are
> > implmented for single system (i mena single processor). So that i
> can
> > parallelize them for multiple processors (distributed as well as
> shared).
> >
> > Thanking in advance for your help.
> >
> > Regards
> > Vaijanath N. Rao
> >
> > _______________________________________________
> > Insight-users mailing list
> > Insight-users@itk.org
> > http://www.itk.org/mailman/listinfo/insight-users
> >
>
>
>
>
> --__--__--
>
> Message: 6
> Date: Tue, 11 May 2004 18:21:34 -0400
> From: Luis Ibanez <luis.ibanez@kitware.com>
> To: Bing Jian <bjian@cise.ufl.edu>
> Cc: Insight Mailing List <Insight-users@itk.org>
> Subject: Re: [Insight-users] Segmentation Programs in ITKApplication
>
>
> Hi Bing,
>
> There are *many* segmentation examples in InsightApplications.
>
> If you are so kind to give us more details regarding the programs
> you are having trouble with, we will be happy to help you get
> around those problems.
>
>
> As you may understand,
> We will need basic information such as:
>
>
> A) The NAME of the program :-)
>
> B) The type of image that you are trying to load.
>
> C) When in the process it is crashing.
>
>
> Please share this information with us,
>
>
> Thanks
>
>
> Luis
>
>
> ---------------------
> Bing Jian wrote:
>
> > Hi, Everyone,
> >
> > I don't know why I cannot run those demo segementation
> > programs with GUI in Insight Application. The first step
> > of loading input image is fine, but the second step will
> > give abort information. That's true for both linux and windows.
> > Anybody can explain it to me? Many thanks in advance!
> >
> >
>
>
>
>
> --__--__--
>
> Message: 7
> Date: Tue, 11 May 2004 19:09:51 -0400 (EDT)
> From: Bing Jian <bjian@cise.ufl.edu>
> To: Luis Ibanez <luis.ibanez@kitware.com>
> Cc: Insight Mailing List <Insight-users@itk.org>
> Subject: Re: [Insight-users] Segmentation Programs in ITKApplication
>
> Hi Luis,
>
> Sorry for not describing my problem clearly.
> Actually I have such problem with almost all segmentation
> GUI demos like GeodesicActiveContour, LevelSetSegmentation,
> FastMarchingLevelSet, etc. The images I am trying to load
> are just 2D *.png images from ITK datasets, I can view them
> after loading. But if I continue the processing, for example,
> click the 'Gradient Magnitude' button in Geodesic Active Contour
> program, it crashes with "Aborted" message. I do not change
> any default parameters.
>
> --
> Best wishes,
> Bing Jian
> bjian@cise.ufl.edu
>
>
> On Tue, 11 May 2004, Luis Ibanez wrote:
>
> >
> > Hi Bing,
> >
> > There are *many* segmentation examples in InsightApplications.
> >
> > If you are so kind to give us more details regarding the programs
> > you are having trouble with, we will be happy to help you get
> > around those problems.
> >
> >
> > As you may understand,
> > We will need basic information such as:
> >
> >
> > A) The NAME of the program :-)
> >
> > B) The type of image that you are trying to load.
> >
> > C) When in the process it is crashing.
> >
> >
> > Please share this information with us,
> >
> >
> > Thanks
> >
> >
> > Luis
> >
> >
> > ---------------------
> > Bing Jian wrote:
> >
> > > Hi, Everyone,
> > >
> > > I don't know why I cannot run those demo segementation
> > > programs with GUI in Insight Application. The first step
> > > of loading input image is fine, but the second step will
> > > give abort information. That's true for both linux and windows.
> > > Anybody can explain it to me? Many thanks in advance!
> > >
> > >
> >
> >
> >
> >
>
>
> --__--__--
>
> Message: 8
> Date: Tue, 11 May 2004 19:13:26 -0400 (EDT)
> From: Bing Jian <bjian@cise.ufl.edu>
> To: Kai Li <likai@cs.uoregon.edu>
> Cc: Insight Mailing List <Insight-users@itk.org>
> Subject: Re: [Insight-users] Segmentation Programs in ITKApplication
>
> Hi, Kai,
>
> Thanks for your information. But I will get "problems reading
> file
> format" when loading a 3D mhd image.
>
>
> --
> Best wishes,
> Bing Jian
> bjian@cise.ufl.edu
>
>
> On Tue, 11 May 2004, Kai Li wrote:
>
> > Hi,
> > I tried the GeodesicActiveContour segmentation application.
> > I guess I may got the similar problems to yours. What I found out
> was that
> > the program works fine for 3D images, but not for 2D images. The
> image
> > reader filter will generate a 3D image (but the size of one of the
> > dimension is 1) even if the input file is a 2D image. I guess it
> is the
> > 3D image that has one dimensiono of size 1 that causes the
> problem.
> >
> > kai
> >
> > On Tue, 11 May 2004, Bing Jian wrote:
> >
> > > Hi, Everyone,
> > >
> > > I don't know why I cannot run those demo segementation
> > > programs with GUI in Insight Application. The first step
> > > of loading input image is fine, but the second step will
> > > give abort information. That's true for both linux and windows.
> > > Anybody can explain it to me? Many thanks in advance!
> > >
> > >
> > > --
> > > Best wishes,
> > > Bing Jian
> > > bjian@cise.ufl.edu
> > >
> > >
> > > _______________________________________________
> > > Insight-users mailing list
> > > Insight-users@itk.org
> > > http://www.itk.org/mailman/listinfo/insight-users
> > >
> > _______________________________________________
> > Insight-users mailing list
> > Insight-users@itk.org
> > http://www.itk.org/mailman/listinfo/insight-users
> >
> >
>
>
>
> --__--__--
>
> _______________________________________________
> Insight-users mailing list
> Insight-users@itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
>
> End of Insight-users Digest
>
More information about the Insight-users
mailing list