Difference between revisions of "ITK Registration Optimization"

From NAMIC Wiki
Jump to: navigation, search
m (Text replacement - "[http://www.na-mic.org/Wiki/images/" to "[https://na-mic.org/w/images/")
 
(75 intermediate revisions by 6 users not shown)
Line 1: Line 1:
−
= Quick Links =
+
{|border=1
 +
|This page contains a record of the development of RegisterImages.  For user documentation please see [https://www.slicer.org/wiki/Modules:RegisterImages-Documentation-3.4 RegisterImages-Documentation-3.4]
 +
|}
  
−
# [http://public.kitware.com/dashboard.php?name=BWHITKOptimization Dashboard for this project]
+
= Slicer3 Module: RegisterImages =
−
# [http://public.kitware.com/dashboard.php?name=BatchMake Dashboard for BatchMake]
 
−
# [http://www.insight-journal.org/batchmake Batchboard (nightly experiment results) for this project]
 
  
−
= Results and Publications =
+
The RegisterImages module is the product of the research discussed below.
  
−
# [http://insight-journal.org/InsightJournalManager/view_reviews.php?pubid=172 Aylward, Stephen; Jomier, Julien; Barre, Sebastien; Davis, Brad; Ibanez, Luis, "Optimizing ITK’s Registration Methods for Multi-processor, Shared-Memory Systems." MICCAI Open Source and Open Data Workshop, 2007] [http://insight-journal.org/InsightJournalManager/download_publication.php?pubid=172&revision=2&name=OptimizingITKRegistrationMethods.pdf&pdf=1 (Download PDF)]
+
== Major Features ==
 +
 
 +
The major features of the module include:
 +
* Default parameters register many full-head and skull-stripped MRI: rigid, affine, and BSpline
 +
* Offers a complete, pipeline-based registration solution
 +
** Load and apply existing transforms
 +
** Compute rigid, affine, and bspline transforms in sequence with a single command
 +
* Intuitive parameters
 +
** Instead of setting obscure "scales" for parameters, you set global values for "Expected Offset", "Expected Rotation", ... to indicate how much mis-registration is anticipated in the data being registered
 +
* MinimizeMemory option provides a way to compute bspline registrations using a dense set of control points and a large number of samples on "normal" computers (albeit computation time increases)
 +
* SampleFromOverlap option allows images of vastly different sizes to be registered
 +
** Helps to avoid (but does not completely eliminate) the annoying ITK exception, "too many samples falls outside of the image"
 +
* Incorporates testing
 +
** Specify a baseline image, and modules will perform the requested registration, compare its results with the baseline image, and return success/failure
 +
* Based on an extensible and re-usable class structure.
 +
 
 +
Each of these features is discussed next.
 +
 
 +
== Head MRI Registration ==
 +
 
 +
=== Example 1: Problem Cases ===
 +
 
 +
* Registers images which were not well resolved (or produced seg-faults) using other slicer registration modules:
 +
** https://www.slicer.org/wiki/Slicer3:Registration
 +
 
 +
==== Difficult affine registration, Same subject, Different protocols: T2 and Fractional Anisotropy ====
 +
{|
 +
|-
 +
| [[Image:RegisterImages-T2-fa-registration-setup-2008-07-30.png|thumb|250px|Setup: T2 (fixed) with FA (moving)]]
 +
| [[Image:RegisterImages-T2-fa-registration-results-2008-07-30.png|thumb|250px|Affine Results: RegisterImages, Default Options (2008-07-30)]]
 +
|-
 +
|}
 +
 
 +
==== Difficult affine registration, Same subject, Different protocols: T1 and Gradient ====
 +
 
 +
{|
 +
|-
 +
| [[Image:RegisterImages-Debrecen-affine-setup-2008-07-30.png|thumb|250px|Setup: T1 (fixed) with Gradient (moving)]]
 +
| [[Image:RegisterImages-Debrecen-affine-results-2008-07-30.png|thumb|250px|Affine Results: RegisterImages, Default Options (2008-07-30)]]
 +
|-
 +
|}
 +
 
 +
=== Example 2: Affine Registration ===
 +
 
 +
* Task:
 +
** Affine registration of head MRI from two different subjects
 +
* Data:
 +
** Using cases UNC-Healthy-Normal002 (fixed) and UNC-Healthy-Normal004 (moving)
 +
** Data provided by Dr. Bullitt at UNC.
 +
** Data is available from Kitware's MIDAS archive at  http://hdl.handle.net/1926/542
 +
** Data can be automatically downloaded into ${RegisterImages_BINARY_DIR}/Testing/Data directory by enabling the CMake variable "BUILD_REGISTER_IMAGES_REAL_WORLD_TESTING"
 +
*** Warning this also enables additional tests that can take 4+ hours to complete.
 +
*** To see the code for automatically downloading from MIDAS (via svn), see Slicer3/Applications/CLI/RegisterImagesModule/Applications/CMakeLists.txt
 +
 
 +
==== Affine registration of head MRI from two difference subjects ====
 +
 
 +
{|
 +
|-
 +
| [[Image:RegisterImages-Normal002-004-affine-setup-2008-07-30.png|thumb|250px|Setup: UNC Normal 002 (fixed) with UNC Normal 004 (moving)]]
 +
| [[Image:RegisterImages-Normal002-004-affine-results-2008-07-30.png|thumb|250px|PipelineAffine Results: RegisterImages, Default Options (2008-07-30)]]
 +
|-
 +
|}
 +
 
 +
==== Affine registration of skull-stripped head MRI from two difference subjects ====
 +
 
 +
{|
 +
|-
 +
| [[Image:RegisterImages-Normal002-004-stripped-affine-setup-2008-07-30.png|thumb|250px|Setup: Skull Stripped: UNC Normal 002 (fixed) with UNC Normal 004 (moving)]]
 +
| [[Image:RegisterImages-Normal002-004-stripped affine-results-2008-07-30.png|thumb|250px|PipelineAffine Results: Skull Stripped RegisterImages, Default Options (2008-07-30)]]
 +
|-
 +
|}
 +
 
 +
=== Example 3: BSpline Registration ===
 +
 
 +
* Task:
 +
** BSpline registration of head MRI from two different subjects
 +
* Data:
 +
** Using cases UNC-Healthy-Normal002 (fixed) and UNC-Healthy-Normal004 (moving)
 +
** Data provided by Dr. Bullitt at UNC.
 +
** Data is available from Kitware's MIDAS archive at  http://hdl.handle.net/1926/542
 +
** Data can be automatically downloaded into ${RegisterImages_BINARY_DIR}/Testing/Data directory by enabling the CMake variable "BUILD_REGISTER_IMAGES_REAL_WORLD_TESTING"
 +
*** Warning this also enables additional tests that can take 4+ hours to complete.
 +
*** To see the code for automatically downloading from MIDAS (via svn), see Slicer3/Applications/CLI/RegisterImagesModule/Applications/CMakeLists.txt
 +
 
 +
==== BSpline, Default Options, Different Subjects, Same Protocol, Head MRI ====
 +
{|
 +
|-
 +
| [[Image:RegisterImages-Normal002-004-bspline-setup-2008-07-30.png|thumb|250px|Setup: UNC Normal 002 (fixed) with UNC Normal 004 (moving)]]
 +
| [[Image:RegisterImages-Normal002-004-bspline-results-2008-07-30.png|thumb|250px|PipelineBSpline Results: RegisterImages, Default Options (2008-07-30)]]
 +
|-
 +
|}
 +
 
 +
==== BSpline, Dense Control-Point Grid, Different Subjects, Same Protocol, Head MRI ====
 +
{|
 +
|-
 +
| [[Image:RegisterImages-Normal002-004-bspline2-setup-2008-07-30.png|thumb|250px|PipelineBSpline Results: UNC Normal 002 (fixed) with UNC Normal 004 (moving), Default Options]]
 +
| [[Image:RegisterImages-Normal002-004-bspline2-results-2008-07-30.png|thumb|250px|PipelineBSpline Results: RegisterImages, Options: MinimizeMemory and ControlPointSpacing=15 pixels (vs. default ControlPointSpacing=40 pixels) (2008-07-30)]]
 +
|-
 +
|}
 +
 
 +
==== BSpline, (Default Options vs. Dense Control-Point Grid), Different Subjects, Same Protocol, Skull-Stripped Head MRI ====
 +
{|
 +
|-
 +
| [[Image:RegisterImages-Normal002-004-stripped-bspline-setup-2008-07-30.png|thumb|250px|Setup: Skull Stripped: UNC Normal 002 (fixed) with UNC Normal 004 (moving)]]
 +
| [[Image:RegisterImages-Normal002-004-stripped-bspline-results-2008-07-30.png|thumb|250px|PipelineBSpline Results: RegisterImages, Default Options (2008-07-30)]]
 +
| [[Image:RegisterImages-Normal002-004-stripped bspline-results2-2008-07-30.png|thumb|250px|PipelineBSpline Results: RegisterImages, MinimizeMemeory and ControlPointSpacing=15 pixels (vs. default ControlPointSpacing=40 pixels) (2008-07-30)]]
 +
|-
 +
|}
 +
 
 +
== Pipeline Registration ==
 +
 
 +
The module implements a registration pipeline.  The steps in that pipeline are as follows:
 +
 
 +
* Step 1: Loaded transform
 +
** You may load a pre-computed transform to initialize the registration.
 +
** If one is loaded, it is immediately applied (i.e., the moving image is resampled)
 +
* Step 2: Initial registration
 +
** Options are:
 +
*** None (sets the center of rotation to the center of the moving image)
 +
*** Landmark (uses N-pairs of landmarks (passed as vectors) and a least-squared error metric to register the images using a rigid transform
 +
*** Image Centers (shifts the images to align their centers)
 +
*** Centers of Mass (shifts the images to align their centers of mass)
 +
*** Second Moments (shifts and rotates the images to align the 1st and 2nd moments)
 +
* Step 3: Registration
 +
** Options are:
 +
*** None (applies the loaded transforms)
 +
*** Initial
 +
**** computes and applies the initial transform to the loaded registrations)
 +
*** Rigid
 +
**** computes a rigid transform and then applies it to the loaded registrations
 +
*** Affine
 +
**** computes an affine transform and then applies it to the loaded registrations
 +
*** BSpline
 +
**** computes a bspline transform and then applies it to the loaded registrations
 +
*** PipelineRigid
 +
**** computes a rigid transform (initialized using the results from the initial registration) and then applies it to the loaded registrations
 +
*** PipelineAffine
 +
**** computes a rigid transform (initialized using the results from the initial registration), uses those results to initialize and compute an affine transform, and then applies it to the loaded registrations
 +
*** PipelineBSpline
 +
**** computes a rigid transform (initialized using the results from the initial registration), uses those results to initialize and compute an affine transform, and then applies it to the loaded registrations, THEN computes and applies a BSpline transform
 +
 
 +
== Intuitive Parameters ==
 +
 
 +
* In rare cases (given unusual acquisition conditions and/or highly inconsistent acquisition protocols) you will need to change the default parameters.
 +
* More often you may wish to tweak parameters to achieve your application-specific speed-vs-accuracy tradeoff
 +
 
 +
=== IO Tab ===
 +
* Set the fixed and moving images using images in the scene
 +
* Optionally set the ResampleImage to store the output image
 +
** If not set, registration won't conduct the final resampling, saving computation time
 +
 
 +
=== Registration Parameters Tab ===
 +
* Load Transform
 +
** provide the Loaded Transform for the loaded phase of registration
 +
* Save Transform
 +
** results of the entire registration pipeline will be saved here
 +
* Initialization
 +
** see registration pipeline discussion
 +
* Registration
 +
** see registration pipeline discussion
 +
** For rigid and affine registrations, one-plus-one evoluation optimization is first applied for N iterations, and then FRPR gradient-line-search optimization is applied.
 +
*** For more information, check the code: RegisterImagesModule/itkOptimizedImageToImageRegistrationmethod.h/txx
 +
** For BSpline registration, a hierarchical registration scheme is used. An image pyramid having 3 levels is used to resample the images and the control grids.  Heuristics are used to control the various resampling parameters.  At each level, registration is conducted using FRPR gradient-line-search optimization.
 +
*** For more information, check the code: RegisterImagesModule/itkBSplineImageToImageRegistrationMethod.h/txx
 +
* Metric
 +
** Use the Mutual Information metric.  It is an multithreaded and optimized version of the Mattes MI method.
 +
*** For more information, check the code; Insight/Code/Review/itkOptMattesMutualInformationImageMetric.h/txx
 +
* "Expected" values
 +
**  For rigid, affine, and bspline registration, parameter scales (refer to the Insight Software Guide) are represented as hyper-parameters in the RegisterImages module.
 +
*** "Expected Offset" controls the offset scales in rigid and affine registration the deformation vector scale in bspline registration
 +
*** "Expected Rotation" is roughly in terms of radians.  It controls the rotation angles in rigid and affine registration
 +
*** "Expected Scale" is for scaling during affine registration
 +
*** "Expected Skew" is for skew for affine registration
 +
 
 +
=== Advaned Registration Parameters Tab ===
 +
* Verbosity level
 +
** Controls the level of detail in the reports in the log file
 +
* Sample from fixed/moving overlap
 +
** When the fixed image is much larger than the moving image, it is CRITICAL to set this flag and to pick a good initialization method.  In that way, only the portion of the fixed image that is initially covered by the moving image will be used during registration.  This prevents ITK from throwing an exception (error) stating that too many fixed-image samples miss (map outside of) the moving image.
 +
* Fixed image intensity percentage threshold
 +
** A less robust way to overcome the image overlap issue discussed above, you can specify a threshold as a portion (0 to 1) of the fixed image intensity range that should be used to select fixed image samples for computing the metric.  That is, by specifying 0.5, only the pixels in the upper half of the fixed-image's intensity range will be used during random sample selection.
 +
** Remember, it is important to include pixels inside and outside of the object of interest, otherwise the fixed image histogram may be too homogeneous for mis-registrations to be detected.
 +
* Random number seed
 +
** To ensure consistent performance, you can set a seed - repeated runs should produce identical results.
 +
* Number of threads
 +
** Number of multi-core/mult-processor threads to use during metric value computations.
 +
* MimimizeMemory
 +
** Turns off caching of intermediate values during bspline registration
 +
** Provides a way to compute bspline registrations using a dense set of control points and a large number of samples on "normal" computers (albeit computation time increases)
 +
** Rule of thumb, if the BSpline registration crashes - re-run with this option enabled.
 +
* use windowed sinc for final interpolation
 +
** If you have time to kill.  Extremely slow and only marginally better than bspline resampling (the default).
 +
 
 +
=== Registration Testing Parameters ===
 +
* Baseline Image
 +
** Set the image against which the Resampled Image (IO tab) will be compared after registration
 +
* Number of Failed Pixels Tolerance
 +
** Registration returns "failure" if this many pixels are different between the Resampled and Baseline images
 +
* Intensity Tolerance
 +
** Minimum intensity difference between corresponding Resampled and Baseline pixels for those pixels to be counted as failures
 +
* Radius Tolerance
 +
** The program will search this neighborhood size about each Resampled pixel to find the closest matching Baseline pixel.   The closest matching pixels are compared using the Intensity Tolerance (above)
 +
* Baseline Difference Image
 +
** Result of subtracting the resampled image from the baseline image
 +
* Baseline Resamples Moving Image
 +
** resampled image, resampled into the space of the baseline image
 +
 
 +
=== Advanced Initial Registration Parameters ===
 +
* Fixed / Moving Landmarks
 +
** A vector string (comma separated base-3 list) of the indexes of corresponding points in the fixed and moving images
 +
** If supplied, then choose "Landmarks" as the initial registration method (see discussion on registration pipeline)
 +
 
 +
=== Advanaced Rigid and Affine Parameters ===
 +
* MaxIterations
 +
** Number of iterations for one-plus-one and for FRPR registration
 +
* Sampling Ratio
 +
** Portion of the image pixels to be used when computing the metric
 +
 
 +
=== Advanced BSpline Parameters ===
 +
* MaxIterations
 +
** Number of iterations for one-plus-one and for FRPR registration
 +
* Sampling Ratio
 +
** Portion of the image pixels to be used when computing the metric
 +
** Do the math...if you have 40 pixels between control points, then there will be 40^3 (64,000) pixels relevant to each control point.  That excessive for directing one control point. Keep the sampling small.  For 40 pixels between control points, a sampling density of 0.1 provide 6,400 pixels for metric computation at each control point - more than enough.
 +
** When in doubt, turn on MinimizeMemory
 +
* Control point spacing (pixels)
 +
** Don't think about grid size - instead think about the level of detail that needs to be resolved (see discussion on sampling ratio).
 +
** When in doubt, turn on MinimizeMemory
 +
 
 +
== Incorporates testing ==
 +
 
 +
* See discussion on the "Registration Testing Parameters" tab.
 +
 
 +
 
 +
 
 +
== Class structure ==
 +
 
 +
* Try it, you'll like it.
 +
* Follows the coding style of itk
 +
* Limited comments, but meaningful variable names
 +
* No documentation is provided or planned - don't even ask.
 +
 
 +
= Instructions for Enabling the RegisterImages module =
 +
 
 +
This module should now be built and distributed by default.  No special steps are needed to use this module.
  
−
* One remaining, high priority task is to complete the integration of the new, threaded, registration methods into ITK.  Luis and Sebastien have adapted the new methods to be 100% backward compatible with ITK's existing classes.  This is a major effort involving approximately 50,000 lines of new code and over 400 new tests in ITK.  The new registration framework is going to be significantly better tested as well as significantly faster than the existing ITK registration framework.  Once it is ported, helper-classes will be added to ITK, and modules using those helper classes will be distributed with Slicer.  We have chosen to spend the time to integrate with ITK because it will serve the broader community, it will benefit from the support of the broader community, it will avoid having to incorporate another SVN checkout into Slicer's build process, and it will keep us from having to maintain and monitor separate dashboards for this effort.
+
= Background =
  
−
= Goals =
+
== Goals ==
  
 
There are two components to this research
 
There are two components to this research
 
# Identify registration algorithms that are suitable for non-rigid registration problems that are endemic to NA-MIC
 
# Identify registration algorithms that are suitable for non-rigid registration problems that are endemic to NA-MIC
−
# Develop implementations of those algorithms that take advantage of multi-core and multi-processor hardware.
+
# Develop implementations of those algorithms that take advantage of multi-core and multi-processor hardware
 +
 
 +
== Steps involved ==
 +
 
 +
# Modify ITK's registration framework to support oriented images
 +
# Modify ITK's registration framework to be thread safe
 +
# Develop multi-threaded versions of select registration modules
 +
# Make everything backward compatible with ITK's existing registration methods and framework
 +
# Deliver in ITK
 +
# Develop helper classes and write IJ article
 +
 
 +
Target date for these deliverables: Jan 1, 2008
 +
 
 +
== Planned follow-on work ==
 +
 
 +
Devise a new metric for MI registration
 +
# If we always use every voxel for the metric, then we can cache the weights by the voxel's position wrt the adjacent control points. For example, for Kilian's situation of a control point every 2 voxels, then there really are only a few unique weight sets that are repeated throughout the volume. Luis had already brought up a variation on this idea.
 +
# This method could also be combined with the rule to not evaluate voxels or control points that fall on background voxels.  This too has been discussed, but such a rule makes multi-threading tricky in that we don't want to waste threads by allocating them to image regions that contain only background voxels.
 +
# The metric could be closely tied to a multiresolution registration scheme.  In fact, the grid and the image resolutions should perhaps be linearly related.  That is, we could tie the metric computation to the resolution of the deformation grid by subsampling the image.  There are situations where this is not a right thing to do (just because the grid is coarse doesn't mean that a small movement isn't important); HOWEVER, as part of a multiresolution registration strategy, it is perhaps the viable option. This would need to be evaluated on the data.
 +
# Have "don't-care" regions in which bspline control points are processed/don't move, e.g., no need to adjust ones that only contain background
 +
 
 +
== Status and News ==
 +
Thanks (but not your questions or comments) go to
 +
* Luis Ibanez, Matt Turek, Stephen Aylward
 +
Questions and comments should go to the Slicer Developers' list
 +
 
 +
== Publications ==
 +
 
 +
# [http://insight-journal.org/InsightJournalManager/view_reviews.php?pubid=172 Aylward, Stephen; Jomier, Julien; Barre, Sebastien; Davis, Brad; Ibanez, Luis, "Optimizing ITK’s Registration Methods for Multi-processor, Shared-Memory Systems." MICCAI Open Source and Open Data Workshop, 2007] [http://insight-journal.org/InsightJournalManager/download_publication.php?pubid=172&revision=2&name=OptimizingITKRegistrationMethods.pdf&pdf=1 (Download PDF)]
 +
# [[NAC_Grid_Enabled_ITK | BWH Neuroimaging Analysis Center (NAC), 2007-2008: Grid Enabled ITK]]
 +
# IJ article on oriented images and registration in ITK
 +
#* http://www.insight-journal.org/dspace/bitstream/1926/1293/2/Brooks_Arbel_FastOrientedImage_V1.pdf
 +
#* Solution presented by the authors is closely related to the changes  made in ITK
  
 
== Algorithmic Requirements and Use Cases ==
 
== Algorithmic Requirements and Use Cases ==
Line 51: Line 327:
 
*# 16 core, Sun Fire, AMDOpteron (UNC: Styner)
 
*# 16 core, Sun Fire, AMDOpteron (UNC: Styner)
  
−
== Data ==
+
== Historic Results ==
−
 
 
−
* Now distributed with CVS
 
−
 
 
−
= Workplan =
 
−
 
 
−
== Establish testing and reporting infrastructure ==
 
−
# Identify timing tools
 
−
## Cross platform and multi-threaded
 
−
## Timing and profiling
 
−
# Develop performance dashboard for collecting results
 
−
## Each test will report time and accuracy to a central server
 
−
## The performance of a test, over time, for a given platform can be viewed on one page
 
−
## The performance of a set of tests, at one point in time, for all platforms can be viewed on one page
 
−
 
 
−
 
 
−
== Develop tests ==
 
−
# Develop modular tests
 
−
# Develop complete registration solutions for use cases
 
−
 
 
−
 
 
−
== ITK Optimization ==
 
−
* Target bottlenecks
 
−
** Multi-thread metric calculation
 
−
*** Initial target is MattesMutualInformationImageToImageMetric
 
−
** Optimize code
 
−
*** Sacrifice some memory and algorithm initialization speed to gain algorithm operation speed increases
 
−
*** Call multi-threaded functions when possible
 
−
* Integrate metrics with transforms and interpolators for tailored performance
 
−
 
 
−
=== MattesMutualInformationImageToImageMetric ===
 
−
==== Optimizations Employed ====
 
−
* GetValue
 
−
** Added multi-threading to GetValue function
 
−
*** Partitions the samples - thereby distributes the computation of the transforms and interpolations across threads
 
−
*** Added the pre-computation of the FixedImageMarginalPDF for the sample to reduce the need for the thread mutex lock
 
−
**** Required the concept of an AdjustedFixedImageMarginalPDF that is updated when a fixed image voxel does not map into the moving image and thereby isn't valid for the current computations.  By only updating when samples are missed, mutex lock to update a cross-thread data structure is needed less often.
 
−
*** Each thread now has its own copy of the joinPDF.  After threads complete, jointPDFs from each thread are summed.  This eliminates mutex from the main loop over samples.
 
−
*** SUMMARY: Speedup on a dual-core system is about 30% (reduction in computation time) when using linear transform and linear interpolation and about 45% when using bspline transform and bspline interpolation.
 
−
** Algorithm optimization (aside from adding multi-threading)
 
−
*** None at this level.  Will be done for transforms and interpolators.
 
−
* GetDerivative
 
−
** Following the same convention as used with GetValue
 
−
* See Publications and Results section for more details
 
−
 
 
−
= Modular tests =
 
−
 
 
−
All tests send two values to performance dashboards
 
−
* the time required
 
−
* an measure of the error (0 = no error; 1 = 100% error)
 
−
 
 
−
Tests being developed and their parameter spaces
 
−
# NearestNeighborInterpTest <numThreads> <dimSize> <factor> [<outputImage>]
 
−
# LinearInterpTest <numThreads> <dimSize> <factor> [<outputImage>]
 
−
# BSplineInterpTest <numThreads> <dimSize> <factor> <bSplineOrder> [<outputImage>]
 
−
# SincInterpTest <numThreads> <dimSize> <factor> [<outputImage>]
 
−
#* Uses the Welch window function
 
−
# BSplineTransformLinearInterpTest <numThreads> <dimSize> <numNodesPerDim> <bSplineOrder> [<outputImage>]
 
−
#* 3 nodes are also added outside of the image for interpolation
 
−
# MeanSquaresImageToImageMetricTest <numThreads> <dimSize> <iterations>
 
−
# CorrelationCoefficientHistogramMetricTest <numThreads> <dimSize> <iterations>
 
−
# NormalizedCorreltationImageToImageMetricTest <numThreads> <dimSize> <iterations>
 
−
# MattesMutualInformationImageToImageMetricTest <numThreads> <dimSize> <numSamples> <iterations>
 
−
#* MattesMutualInformationMetric defaults to BSpline interpolator - this test overrides to use nearest neighbor interpolation
 
−
# MutualInformationImageToImageMetricTest <numThreads> <dimSize> <numSamples> <iterations>
 
−
# NormalizedMutualInformationHistrogramMetricTest <numThreads> <dimSize> <iterations>
 
−
 
 
−
SECOND GENERATION TEST
 
−
* Computes runtime for GetValue, GetDerivative, and GetValueAndDerivative for standard ITK implementation and for the optimized version being developed.  Also computes difference (if any) between their answers and reports as an error measure
 
−
 
 
−
= Larger tests =
 
  
−
These are top down tests that implement registration tasks as a user would.  These tests are run on real medical image data.  The goal is to use these larger, realistic tests to focus our optimization effort.
+
January 5, 2008 - Note: "Opt" results are not using the OptLinearInterpolateImageFunction.
  
−
*AlignMomentsTest <fixed image> <moving image> [number of threads] [transformed moving image] [post-registration difference image] [pre-registration difference image]
+
*[https://na-mic.org/w/images/f/fe/MattesGetValue.pdf MattesMI GetValue Results]
−
*AlignRigidMSETest <fixed image> <moving image> [number of threads] [number of iterations] [transformed moving image] [post-registration difference image] [pre-registration difference image]
+
*[https://na-mic.org/w/images/1/1b/MattesBSplineGetValue.pdf MattesMI, b-spline interpolation and transform, GetValue Results]
−
*AlignRigidMMITest <fixed image> <moving image> [number of threads] [number of iterations] [transformed moving image] [post-registration difference image] [pre-registration difference image]
+
*[https://na-mic.org/w/images/1/16/MeanSquaresGetValue.pdf MeanSquares GetValue Results]
 +
*[http://wiki.na-mic.org/Wiki/images/c/c5/MattesGetValueAndDerivative.pdf MattesMI GetValueAndDerivative Results]
 +
*[http://wiki.na-mic.org/Wiki/images/f/f6/MattesBSplineGetValueAndDerivative.pdf MattesMI, b-spline interpolation and transform, GetValueAndDerivative Results]
 +
*[http://wiki.na-mic.org/Wiki/images/2/2d/MeanSquaresGetValueAndDerivative.pdf MeanSquares GetValueAndDerivative Results]
  
−
= Events =
+
== Historic Events ==
  
 
* [http://www.na-mic.org/Wiki/index.php/ITK_Registration_Optimization/2007-04-06-tcon April 6, 2007: TCon]
 
* [http://www.na-mic.org/Wiki/index.php/ITK_Registration_Optimization/2007-04-06-tcon April 6, 2007: TCon]
Line 136: Line 345:
 
* [http://www.na-mic.org/Wiki/index.php/ITK_Registration_Optimization/2007-05-01-tcon May 1, 2007: TCon]
 
* [http://www.na-mic.org/Wiki/index.php/ITK_Registration_Optimization/2007-05-01-tcon May 1, 2007: TCon]
 
* [http://www.na-mic.org/Wiki/index.php/ITK_Registration_Optimization/2007-06-27-tcon June 27, 2007: NAMIC Programmers' Week]
 
* [http://www.na-mic.org/Wiki/index.php/ITK_Registration_Optimization/2007-06-27-tcon June 27, 2007: NAMIC Programmers' Week]
 +
* [http://wiki.na-mic.org/Wiki/index.php/Registration_Update January, 2008: NAMIC AHM]
  
 
= Related Pages =
 
= Related Pages =
  
 
* [[Non Rigid Registration]]
 
* [[Non Rigid Registration]]
−
* [[Slicer3:Performance_Analysis]]
+
* [https://www.slicer.org/wiki/Slicer3:Performance_Analysis Slicer3:Performance_Analysis]
 
* [http://www.na-mic.org/Wiki/index.php/User:Barre/ITK_Registration_Optimization User:Barre/ITK Registration Optimization]
 
* [http://www.na-mic.org/Wiki/index.php/User:Barre/ITK_Registration_Optimization User:Barre/ITK Registration Optimization]
 
* [[ITK_Registration_Optimization/Testing_And_Backward_Forward_Compatibility | Testing and ITK Backward Forward Compatibility]]
 
* [[ITK_Registration_Optimization/Testing_And_Backward_Forward_Compatibility | Testing and ITK Backward Forward Compatibility]]
  
−
= Performance Measurement =
+
= Notes on Software Profiling Tools =
 
* [http://www.lw-tech.com/index.php LTProf - simple profilter for Windows - Shareware]
 
* [http://www.lw-tech.com/index.php LTProf - simple profilter for Windows - Shareware]
 
* [http://www.intel.com/cd/software/products/asmo-na/eng/vtune/vlin/239145.htm Intel's VTune for Linux] ($)
 
* [http://www.intel.com/cd/software/products/asmo-na/eng/vtune/vlin/239145.htm Intel's VTune for Linux] ($)

Latest revision as of 18:27, 10 July 2017

Home < ITK Registration Optimization
This page contains a record of the development of RegisterImages. For user documentation please see RegisterImages-Documentation-3.4

Contents

Slicer3 Module: RegisterImages

The RegisterImages module is the product of the research discussed below.

Major Features

The major features of the module include:

  • Default parameters register many full-head and skull-stripped MRI: rigid, affine, and BSpline
  • Offers a complete, pipeline-based registration solution
    • Load and apply existing transforms
    • Compute rigid, affine, and bspline transforms in sequence with a single command
  • Intuitive parameters
    • Instead of setting obscure "scales" for parameters, you set global values for "Expected Offset", "Expected Rotation", ... to indicate how much mis-registration is anticipated in the data being registered
  • MinimizeMemory option provides a way to compute bspline registrations using a dense set of control points and a large number of samples on "normal" computers (albeit computation time increases)
  • SampleFromOverlap option allows images of vastly different sizes to be registered
    • Helps to avoid (but does not completely eliminate) the annoying ITK exception, "too many samples falls outside of the image"
  • Incorporates testing
    • Specify a baseline image, and modules will perform the requested registration, compare its results with the baseline image, and return success/failure
  • Based on an extensible and re-usable class structure.

Each of these features is discussed next.

Head MRI Registration

Example 1: Problem Cases

Difficult affine registration, Same subject, Different protocols: T2 and Fractional Anisotropy

Setup: T2 (fixed) with FA (moving)
Affine Results: RegisterImages, Default Options (2008-07-30)

Difficult affine registration, Same subject, Different protocols: T1 and Gradient

Setup: T1 (fixed) with Gradient (moving)
Affine Results: RegisterImages, Default Options (2008-07-30)

Example 2: Affine Registration

  • Task:
    • Affine registration of head MRI from two different subjects
  • Data:
    • Using cases UNC-Healthy-Normal002 (fixed) and UNC-Healthy-Normal004 (moving)
    • Data provided by Dr. Bullitt at UNC.
    • Data is available from Kitware's MIDAS archive at http://hdl.handle.net/1926/542
    • Data can be automatically downloaded into ${RegisterImages_BINARY_DIR}/Testing/Data directory by enabling the CMake variable "BUILD_REGISTER_IMAGES_REAL_WORLD_TESTING"
      • Warning this also enables additional tests that can take 4+ hours to complete.
      • To see the code for automatically downloading from MIDAS (via svn), see Slicer3/Applications/CLI/RegisterImagesModule/Applications/CMakeLists.txt

Affine registration of head MRI from two difference subjects

Setup: UNC Normal 002 (fixed) with UNC Normal 004 (moving)
PipelineAffine Results: RegisterImages, Default Options (2008-07-30)

Affine registration of skull-stripped head MRI from two difference subjects

Setup: Skull Stripped: UNC Normal 002 (fixed) with UNC Normal 004 (moving)
PipelineAffine Results: Skull Stripped RegisterImages, Default Options (2008-07-30)

Example 3: BSpline Registration

  • Task:
    • BSpline registration of head MRI from two different subjects
  • Data:
    • Using cases UNC-Healthy-Normal002 (fixed) and UNC-Healthy-Normal004 (moving)
    • Data provided by Dr. Bullitt at UNC.
    • Data is available from Kitware's MIDAS archive at http://hdl.handle.net/1926/542
    • Data can be automatically downloaded into ${RegisterImages_BINARY_DIR}/Testing/Data directory by enabling the CMake variable "BUILD_REGISTER_IMAGES_REAL_WORLD_TESTING"
      • Warning this also enables additional tests that can take 4+ hours to complete.
      • To see the code for automatically downloading from MIDAS (via svn), see Slicer3/Applications/CLI/RegisterImagesModule/Applications/CMakeLists.txt

BSpline, Default Options, Different Subjects, Same Protocol, Head MRI

Setup: UNC Normal 002 (fixed) with UNC Normal 004 (moving)
PipelineBSpline Results: RegisterImages, Default Options (2008-07-30)

BSpline, Dense Control-Point Grid, Different Subjects, Same Protocol, Head MRI

PipelineBSpline Results: UNC Normal 002 (fixed) with UNC Normal 004 (moving), Default Options
PipelineBSpline Results: RegisterImages, Options: MinimizeMemory and ControlPointSpacing=15 pixels (vs. default ControlPointSpacing=40 pixels) (2008-07-30)

BSpline, (Default Options vs. Dense Control-Point Grid), Different Subjects, Same Protocol, Skull-Stripped Head MRI

Setup: Skull Stripped: UNC Normal 002 (fixed) with UNC Normal 004 (moving)
PipelineBSpline Results: RegisterImages, Default Options (2008-07-30)
PipelineBSpline Results: RegisterImages, MinimizeMemeory and ControlPointSpacing=15 pixels (vs. default ControlPointSpacing=40 pixels) (2008-07-30)

Pipeline Registration

The module implements a registration pipeline. The steps in that pipeline are as follows:

  • Step 1: Loaded transform
    • You may load a pre-computed transform to initialize the registration.
    • If one is loaded, it is immediately applied (i.e., the moving image is resampled)
  • Step 2: Initial registration
    • Options are:
      • None (sets the center of rotation to the center of the moving image)
      • Landmark (uses N-pairs of landmarks (passed as vectors) and a least-squared error metric to register the images using a rigid transform
      • Image Centers (shifts the images to align their centers)
      • Centers of Mass (shifts the images to align their centers of mass)
      • Second Moments (shifts and rotates the images to align the 1st and 2nd moments)
  • Step 3: Registration
    • Options are:
      • None (applies the loaded transforms)
      • Initial
        • computes and applies the initial transform to the loaded registrations)
      • Rigid
        • computes a rigid transform and then applies it to the loaded registrations
      • Affine
        • computes an affine transform and then applies it to the loaded registrations
      • BSpline
        • computes a bspline transform and then applies it to the loaded registrations
      • PipelineRigid
        • computes a rigid transform (initialized using the results from the initial registration) and then applies it to the loaded registrations
      • PipelineAffine
        • computes a rigid transform (initialized using the results from the initial registration), uses those results to initialize and compute an affine transform, and then applies it to the loaded registrations
      • PipelineBSpline
        • computes a rigid transform (initialized using the results from the initial registration), uses those results to initialize and compute an affine transform, and then applies it to the loaded registrations, THEN computes and applies a BSpline transform

Intuitive Parameters

  • In rare cases (given unusual acquisition conditions and/or highly inconsistent acquisition protocols) you will need to change the default parameters.
  • More often you may wish to tweak parameters to achieve your application-specific speed-vs-accuracy tradeoff

IO Tab

  • Set the fixed and moving images using images in the scene
  • Optionally set the ResampleImage to store the output image
    • If not set, registration won't conduct the final resampling, saving computation time

Registration Parameters Tab

  • Load Transform
    • provide the Loaded Transform for the loaded phase of registration
  • Save Transform
    • results of the entire registration pipeline will be saved here
  • Initialization
    • see registration pipeline discussion
  • Registration
    • see registration pipeline discussion
    • For rigid and affine registrations, one-plus-one evoluation optimization is first applied for N iterations, and then FRPR gradient-line-search optimization is applied.
      • For more information, check the code: RegisterImagesModule/itkOptimizedImageToImageRegistrationmethod.h/txx
    • For BSpline registration, a hierarchical registration scheme is used. An image pyramid having 3 levels is used to resample the images and the control grids. Heuristics are used to control the various resampling parameters. At each level, registration is conducted using FRPR gradient-line-search optimization.
      • For more information, check the code: RegisterImagesModule/itkBSplineImageToImageRegistrationMethod.h/txx
  • Metric
    • Use the Mutual Information metric. It is an multithreaded and optimized version of the Mattes MI method.
      • For more information, check the code; Insight/Code/Review/itkOptMattesMutualInformationImageMetric.h/txx
  • "Expected" values
    • For rigid, affine, and bspline registration, parameter scales (refer to the Insight Software Guide) are represented as hyper-parameters in the RegisterImages module.
      • "Expected Offset" controls the offset scales in rigid and affine registration the deformation vector scale in bspline registration
      • "Expected Rotation" is roughly in terms of radians. It controls the rotation angles in rigid and affine registration
      • "Expected Scale" is for scaling during affine registration
      • "Expected Skew" is for skew for affine registration

Advaned Registration Parameters Tab

  • Verbosity level
    • Controls the level of detail in the reports in the log file
  • Sample from fixed/moving overlap
    • When the fixed image is much larger than the moving image, it is CRITICAL to set this flag and to pick a good initialization method. In that way, only the portion of the fixed image that is initially covered by the moving image will be used during registration. This prevents ITK from throwing an exception (error) stating that too many fixed-image samples miss (map outside of) the moving image.
  • Fixed image intensity percentage threshold
    • A less robust way to overcome the image overlap issue discussed above, you can specify a threshold as a portion (0 to 1) of the fixed image intensity range that should be used to select fixed image samples for computing the metric. That is, by specifying 0.5, only the pixels in the upper half of the fixed-image's intensity range will be used during random sample selection.
    • Remember, it is important to include pixels inside and outside of the object of interest, otherwise the fixed image histogram may be too homogeneous for mis-registrations to be detected.
  • Random number seed
    • To ensure consistent performance, you can set a seed - repeated runs should produce identical results.
  • Number of threads
    • Number of multi-core/mult-processor threads to use during metric value computations.
  • MimimizeMemory
    • Turns off caching of intermediate values during bspline registration
    • Provides a way to compute bspline registrations using a dense set of control points and a large number of samples on "normal" computers (albeit computation time increases)
    • Rule of thumb, if the BSpline registration crashes - re-run with this option enabled.
  • use windowed sinc for final interpolation
    • If you have time to kill. Extremely slow and only marginally better than bspline resampling (the default).

Registration Testing Parameters

  • Baseline Image
    • Set the image against which the Resampled Image (IO tab) will be compared after registration
  • Number of Failed Pixels Tolerance
    • Registration returns "failure" if this many pixels are different between the Resampled and Baseline images
  • Intensity Tolerance
    • Minimum intensity difference between corresponding Resampled and Baseline pixels for those pixels to be counted as failures
  • Radius Tolerance
    • The program will search this neighborhood size about each Resampled pixel to find the closest matching Baseline pixel. The closest matching pixels are compared using the Intensity Tolerance (above)
  • Baseline Difference Image
    • Result of subtracting the resampled image from the baseline image
  • Baseline Resamples Moving Image
    • resampled image, resampled into the space of the baseline image

Advanced Initial Registration Parameters

  • Fixed / Moving Landmarks
    • A vector string (comma separated base-3 list) of the indexes of corresponding points in the fixed and moving images
    • If supplied, then choose "Landmarks" as the initial registration method (see discussion on registration pipeline)

Advanaced Rigid and Affine Parameters

  • MaxIterations
    • Number of iterations for one-plus-one and for FRPR registration
  • Sampling Ratio
    • Portion of the image pixels to be used when computing the metric

Advanced BSpline Parameters

  • MaxIterations
    • Number of iterations for one-plus-one and for FRPR registration
  • Sampling Ratio
    • Portion of the image pixels to be used when computing the metric
    • Do the math...if you have 40 pixels between control points, then there will be 40^3 (64,000) pixels relevant to each control point. That excessive for directing one control point. Keep the sampling small. For 40 pixels between control points, a sampling density of 0.1 provide 6,400 pixels for metric computation at each control point - more than enough.
    • When in doubt, turn on MinimizeMemory
  • Control point spacing (pixels)
    • Don't think about grid size - instead think about the level of detail that needs to be resolved (see discussion on sampling ratio).
    • When in doubt, turn on MinimizeMemory

Incorporates testing

  • See discussion on the "Registration Testing Parameters" tab.


Class structure

  • Try it, you'll like it.
  • Follows the coding style of itk
  • Limited comments, but meaningful variable names
  • No documentation is provided or planned - don't even ask.

Instructions for Enabling the RegisterImages module

This module should now be built and distributed by default. No special steps are needed to use this module.

Background

Goals

There are two components to this research

  1. Identify registration algorithms that are suitable for non-rigid registration problems that are endemic to NA-MIC
  2. Develop implementations of those algorithms that take advantage of multi-core and multi-processor hardware

Steps involved

  1. Modify ITK's registration framework to support oriented images
  2. Modify ITK's registration framework to be thread safe
  3. Develop multi-threaded versions of select registration modules
  4. Make everything backward compatible with ITK's existing registration methods and framework
  5. Deliver in ITK
  6. Develop helper classes and write IJ article

Target date for these deliverables: Jan 1, 2008

Planned follow-on work

Devise a new metric for MI registration

  1. If we always use every voxel for the metric, then we can cache the weights by the voxel's position wrt the adjacent control points. For example, for Kilian's situation of a control point every 2 voxels, then there really are only a few unique weight sets that are repeated throughout the volume. Luis had already brought up a variation on this idea.
  2. This method could also be combined with the rule to not evaluate voxels or control points that fall on background voxels. This too has been discussed, but such a rule makes multi-threading tricky in that we don't want to waste threads by allocating them to image regions that contain only background voxels.
  3. The metric could be closely tied to a multiresolution registration scheme. In fact, the grid and the image resolutions should perhaps be linearly related. That is, we could tie the metric computation to the resolution of the deformation grid by subsampling the image. There are situations where this is not a right thing to do (just because the grid is coarse doesn't mean that a small movement isn't important); HOWEVER, as part of a multiresolution registration strategy, it is perhaps the viable option. This would need to be evaluated on the data.
  4. Have "don't-care" regions in which bspline control points are processed/don't move, e.g., no need to adjust ones that only contain background

Status and News

Thanks (but not your questions or comments) go to

  • Luis Ibanez, Matt Turek, Stephen Aylward

Questions and comments should go to the Slicer Developers' list

Publications

  1. Aylward, Stephen; Jomier, Julien; Barre, Sebastien; Davis, Brad; Ibanez, Luis, "Optimizing ITK’s Registration Methods for Multi-processor, Shared-Memory Systems." MICCAI Open Source and Open Data Workshop, 2007 (Download PDF)
  2. BWH Neuroimaging Analysis Center (NAC), 2007-2008: Grid Enabled ITK
  3. IJ article on oriented images and registration in ITK

Algorithmic Requirements and Use Cases

  • Requirements
    1. relatively robust, with few parameters to tweak
    2. runs on grey scale images
    3. has already been published
    4. relatively fast (ideally speaking a few minutes for volume to volume).
    5. not patented
    6. can be implemented in ITK and parallelized.

Hardware Platform Requirements and Use Cases

  • Requirements
    1. Shared memory
    2. Single and multi-core machines
    3. Single and multi-processor machines
    4. AMD and Intel - Windows, Linux, and SunOS
  • Use-cases
    1. Intel Core2Duo
    2. Intel quad-core Xeon processors, Visual Studio 8, Windows Vista (Kitware: redwall)
    3. 6 CPU Sun, Solaris 8 (SPL: vision)
    4. 12 CPU Sun, Solaris 8 (SPL: forest and ocean)
    5. 16 core Opteron (SPL: john, ringo, paul, george)
    6. 16 core, Sun Fire, AMDOpteron (UNC: Styner)

Historic Results

January 5, 2008 - Note: "Opt" results are not using the OptLinearInterpolateImageFunction.

Historic Events

Related Pages

Notes on Software Profiling Tools

Retrieved from "https://www.na-mic.org/w/index.php?title=ITK_Registration_Optimization&oldid=97801"