ITK Registration Methods Kitware Inc Overview Image Resampling

ITK Registration Methods Kitware Inc.

Overview Image Resampling Registration Framework Multi-Modality Multi-Resolution Deformable registration

Image Resampling

Why Resampling ? Resampling is the Essence of Intensity-Based Image Registration

What is an Image ? An Image is a sampling of a continuous field using a discrete grid

Image Origin & Spacing (Sx) y x Spacing (Sy) Origin (Ox, Oy)

Image Sampling Grid Spacing (Sx) Spacing (Sy) Origin (Ox, Oy)

Image Pixel Spacing (Sx) Pixel Value Pixel Region Spacing (Sy) Origin (Ox, Oy)

Image Indices Spacing (Sx) [0, 7] Pixel Index [4, 7] [0, 6] [0, 5] [0, 4] [0, 3] [0, 2] Spacing (Sy) [0, 1] Origin (Ox, Oy) [0, 0] [1, 0] [2, 0] [3, 0] [4, 0] [5, 0]

Index to Physical Coordinates Spacing (Sx) [0, 7] Pixel Index [4, 7] [0, 6] P[0] = Index[0] x[0, 5] Spacing[0] + Origin[0] P[1] = Index[1] x[0, 4] Spacing[1] + Origin[1] Index[0] = floor( [0, 3] ( P[0] - Origin[0] ) / Spacing[0] + 0. 5 ) Index[1] = floor( ( P[1] - Origin[1] ) / Spacing[1] + 0. 5 ) [0, 2] Spacing (Sy) [0, 1] Origin (Ox, Oy) [0, 0] [1, 0] [2, 0] [3, 0] [4, 0] [5, 0]

Image Region Spacing (Sx) Pixel Value Image Region Pixel Region Spacing (Sy) Origin (Ox, Oy)

Image Region Spacing (Sx) Image Region Size [3, 5] Starting Index [2, 3] Spacing (Sy) Origin (Ox, Oy)

Basic Resampling Trivial Cases

Sub-Sampling by Half Spacing (Sx) Image Region Spacing ( 2 x Sx ) Spacing ( 2 x Sy ) Spacing (Sy) Origin (Ox, Oy)

Sub-Sampling by Half New Spacing S’x New Origin (O’x, O’y) Origin (Ox, Oy) New Spacing S’y

Super-Sampling by Double Spacing (Sx) Image Region Spacing ( Sx/2 ) Spacing ( Sy/2 ) Spacing (Sy) Origin (Ox, Oy)

Super-Sampling by Double New Spacing S’x New Origin (O’x, O’y) New Spacing Origin (Ox, Oy) S’y

Resampling in ITK itk: : Resample. Image. Filter

Resampling in ITK Origin Spacing Region Size Resample Filter Input Image Transform Region Start Interpolator Output Image

Resample Image Filter #include "itk. Image. h" #include "itk. Resample. Image. Filter. h" #include "itk. Identity. Transform. h“ #include "itk. Linear. Interpolate. Image. Function. h" typedef itk: : Image< char, 2 > Image. Type; Image. Type: : Pointer input. Image = Get. Image. Some. How(); typedef itk: : Resample. Image. Filter< Image. Type > Filter. Type; Filter. Type: : Pointer resampler = Filter. Type: : New(); Image. Type: : Size. Type size; size[0] = 200; size[1] = 300; Image. Type: : Index. Type start; start[0] = 0; start[1] = 0;

Resample Image Filter Image. Type: : Point. Type origin; origin[0] = 10. 0; // millimeters origin[1] = 25. 5; // millimeters Image. Type: : Spacing. Type spacing; spacing[0] = 2. 0; // millimeters spacing[1] = 1. 5; // millimeters resampler->Set. Output. Spacing( spacing ); resampler->Set. Output. Origin( origin ); resampler->Set. Size( size ); resampler->Set. Output. Start. Index( start ); resampler->Set. Default. Pixel. Value( 100 ); resampler->Set. Input( input. Image );

Resample Image Filter typedef itk: : Linear. Interpolate. Image. Function< Image. Type, double > Interpolator. Type; Interpolator. Type: : Pointerpolator = Interpolator. Type: : New(); typedef itk: : Translation. Transform< double, 2 > Transform. Type; Transform. Type: : Pointer transform = Transform. Type: : New(); transform->Set. Identity(); resampler->Set. Interpolator( interpolator ); resampler->Set. Transform( transform ); resampler->Update(); const Image. Type * output. Image = resampler->Get. Output();

Basic Registration

j Coordinate System Conversionsj Image Grid i Image Grid y i y’ y x x Physical Coordinates Space Transform x’ Physical Coordinates

Things I will not do… I will not register images in pixel space I will not register images in pixel space I will not register images in pix

Fixed Image & Moving Image j j Fixed Image Grid i Moving Image Grid y y’ x Fixed Image Physical Coordinates i Space Transform x’ Moving Image Physical Coordinates

Selecting Moving & Fixed Images In principle the denomination of Fixed Image & Moving Image is arbitrary In practice the moving image is the one that will be resampled into the fixed image coordinate system

Quiz #1 Images from the same patient MRI-T 2 PET Moving Image ? Fixed Image ? 256 x 256 pixels 128 x 128 pixels Images provided as part of the project: “Retrospective Image Registration Evaluation”, NIH, Project No. 8 R 01 EB 002124 -03, Principal Investigator, J. Michael Fitzpatrick, Vanderbilt University, Nashville, TN.

Quiz #2 Images from the same patient MRI-T 2 PET Scaling Transform What scale factor ? a) 2. 0 b) 1. 0 256 x 256 pixels c) 0. 5 128 x 128 pixels Images provided as part of the project: “Retrospective Image Registration Evaluation”, NIH, Project No. 8 R 01 EB 002124 -03, Principal Investigator, J. Michael Fitzpatrick, Vanderbilt University, Nashville, TN.

Things I will not do… I will not register images in pixel space I will not register images in pixel space I will not register images in pix

Registration in ITK Multi Resolution Registration Framework Image Registration Framework Components PDE Based Registration FEM Based Registration

Registration Framework

Components Registration Method Fixed Image Metric Interpolator Moving Image Transform Optimizer

Image Metrics Mean Squares Normalized Correlation Mean Reciprocal Square Difference Mutual Information - Viola-Wells - Mattes - Histogram based - Histogram normalized

Transforms Translation Scaling Rotation Rigid 3 D Rigid 2 D Affine BSplines: TPS, EBS, VS

Optimizers Gradient Descent Regular Step Gradient Descent Conjugate Gradient Levenberg-Marquardt One plus One Evolutionary Algorithm

Interpolators Nearest Neighbor Linear BSpline

Exercise 20

Image Registration Fixed Image Moving Image Registered Moving Image

Image Registration #include ”itk. Image. Registration. Method. h” #include ”itk. Translation. Transform. h” #include ”itk. Mean. Squares. Image. To. Image. Metric. h” #include ”itk. Linear. Interpolate. Image. Function. h” #include ”itk. Regular. Step. Gradient. Descent. Optimizer. h” #include ”itk. Image. File. Reader. h” #include ”itk. Image. File. Writer. h” #include ”itk. Resample. Image. Filter. h”

Image Registration const unsigned int Dimension = 2; typedef unsigned char Pixel. Type; typedef typedef itk: : Image< Pixel. Type , Dimension > itk: : Translation. Transform< double, Dimension > itk: : Regular. Step. Gradient. Descent. Optimizer itk: : Linear. Interpolate. Image. Function< Moving. Image. Type , double > typedef itk: : Mean. Squares. Image. To. Image. Metric< Fixed. Image. Type , Moving. Image. Type > typedef itk: : Image. Registration. Method< Fixed. Image. Type , Moving. Image. Type > Fixed. Image. Type; Moving. Image. Type; Transform. Type; Optimizer. Type; Interpolator. Type; Metric. Type; Registration. Type;

Image Registration Transform. Type: : Pointer Optimizer. Type: : Pointer Interpolator. Type: : Pointer Metric. Type: : Pointer Registration. Type: : Pointer transform optimizer interpolator metric registrator->Set. Transform( registrator->Set. Optimizer( registrator->Set. Interpolator( registrator->Set. Metric( = Transform. Type: : New(); = Optimizer. Type: : New(); = Interpolator. Type: : New(); = Metric. Type: : New(); = Registration. Type: : New(); transform ); optimizer ); interpolator ); metric ); registrator->Set. Fixed. Image( fixed. Image. Reader->Get. Output() ); registrator->Set. Moving. Image( moving. Image. Reader->Get. Output() );

Image Registration registrator->Set. Fixed. Image. Region( fixed. Image. Reader->Get. Output()->Get. Buffered. Region() ); typedef Registration. Type: : Parameters. Type; transform->Set. Identity(); registrator->Set. Initial. Transform. Parameters( transform->Get. Parameters() ); optimizer->Set. Maximum. Step. Length( 4. 00 ); optimizer->Set. Minimum. Step. Length( 0. 01 ); optimizer->Set. Number. Of. Iterations( 100 ); optimizer->Maximize. Off();

Image Registration try { registrator->Start. Registration (); } catch( itk: : Exception. Object & excp ) { std: : cerr << “Error in registration” << std: : endl; std: : cerr << excp << std: : endl; } transform->Set. Parameters( registrator->Get. Last. Transform. Parameters() );

Image Registration typedef itk: : Resample. Image. Filter< Fixed. Image. Type , Moving. Image. Type > Resampler. Type; Resampler. Type : : Pointer resampler = Resampler. Type: : New(); resampler->Set. Transform ( transform ); resampler->Set. Input( moving. Image. Reader->Get. Output() ); Fixed. Image. Type : : Pointer fixed. Image = fixed. Image. Reader->Get. Output(); resampler->Set. Origin( fixed. Image->Get. Origin() ); resampler->Set. Spacing( fixed. Image->Get. Spacing() ); resampler->Set. Size( fixed. Image->Get. Largest. Possible. Region()->Get. Size() ); resampler->Update();

Tracking the Registration Process

Observing Registration #include ”itk. Command. h” class Command. Iteration : public itk: : Command { public: typedef Command. Iteration Self; typedef itk: : Command Super. Class; typedef itk: : Smart. Pointer< Self > Pointer; itk. New. Macro( Self ); protected: Command. Iteration() {}; public: typedef itk: : Regular. Step. Gradient. Descent. Optimizer. Type; typedef const Optimizer. Type * Optimizer. Pointer;

Observing Registration void Execute( itk: : Object * caller, const itk: : Event. Object & event ) { this-> Execute( (const itk: : Object *) caller, event ); } void Execute( const itk: : Object * caller, const itk: : Event. Object & event ) { Optimizer. Pointer optimizer = dynamic_cast< Optimizer. Pointer >( caller ); if( typeid( event ) == typeid( itk: : Iteration. Event ) ) { std: : cout << optimizer->Get. Current. Iteration() << “ : “; std: : cout << optimizer->Get. Value() << “ : “; std: : cout << optimizer->Get. Current. Position() << std: : endl; } }

Image Registration Command. Iteration: : Pointer observer = Command. Iteration: : New(); optimizer->Add. Observer( itk: : Iteration. Event(), observer ) try { registrator->Start. Registration (); } catch( itk: : Exception. Object & excp ) { std: : cerr << “Error in registration” << std: : endl; std: : cerr << excp << std: : endl; }

Image Similarity Metrics

Image Metrics How similar is image A to image B ?

Image Metrics Does Image B matches Image A better than Image C ?

Image Metrics Match( A , B ) Image B > < Image A Match( A , C ) Image C

Image Metrics Match( A , B ) Simplest Metric Mean Squared Differences

Mean Squared Differences For each pixel in A Image B Difference( index ) = A( index ) – B( index ) Sum += Difference( index ) 2 Match( A , B ) = Sum / number. Of. Pixels

j For each pixel in the Fixed Image j Fixed Image Grid i Moving Image Grid y y’ x Fixed Image Physical Coordinates i Space Transform x’ Moving Image Physical Coordinates

Image Metrics Fixed Image Valu e Metric Moving Image Interpolator Transform Parameters

Mean Squared Differences #include "itk. Image. h" #include "itk. Mean. Squares. Image. To. Image. Metric. h" #include "itk. Linear. Interpolate. Image. Function. h" #include "itk. Translation. Transform. h" typedef itk: : Image< char, 2 > Image. Type; Image. Type: : Const. Pointer fixed. Image = Get. Fixed. Image(); Image. Type: : Const. Pointer moving. Image = Get. Moving. Image(); typedef itk: : Linear. Interpolate. Image. Function< Image. Type, double > Interpolator. Type; Interpolator. Type: : Pointerpolator = Interpolator. Type: : New(); typedef itk: : Translation. Transform< double, 2 > Transform. Type; Transform. Type: : Pointer transform = Transform. Type: : New();

Mean Squared Differences typedef itk: : Mean. Squares. Image. To. Image. Metric< Image. Type, Image. Type > Metric. Type; Metric. Type: : Pointer metric = Metric. Type: : New(); metric->Set. Interpolator( interpolator ); metric->Set. Transform( transform ); metric->Set. Fixed. Image( fixed. Image ); metric->Set. Moving. Image( moving. Image ); Metric. Type: : Transform. Parameters. Type translation( Dimension ); translation[0] = 12; translation[1] = 27; double value = metric->Get. Value( translation );

Mean Squared Differences Metric. Type: : Transform. Parameters. Type translation( Dimension ); double value[21]; for( int dx = 0; dx <= 20; dx++) { for( int dy = 0; dy <= 20; dy++) { translation[0] = dx; translation[1] = dy; value[dx][dy] = metric->Get. Value( translation ); } }

Evaluating many matches y y Transform x Fixed Image x Moving Image

Plotting the Metric Mean Squared Differences Transform Parametric Space

Plotting the Metric Mean Squared Differences Transform Parametric Space

Evaluating many matches y y Transform (-15, -25) mm x Fixed Image x Moving Image

Plotting the Metric Mean Squared Differences Transform Parametric Space

Plotting the Metric Mean Squared Differences Transform Parametric Space

The Best Transform Parameters Evaluation of the full parameter space is equivalent to performing optimization by exhaustive search

The Best Transform Parameters Very Safe but Very Slow

The Best Transform Parameters Better Optimization Methods for example Gradient Descent

Image Registration Framework Fixed Image Metric Moving Image Interpolator Optimizer Transform Parameters

Gradient Descent Optimizer f( x , y ) S = Step L = Learning Rate ∆ G( x , y ) = f( x , y ) S = L ∙ G( x , y )

Gradient Descent Optimizer f( x , y ) ∆ G( x , y ) = f( x , y ) S = L ∙ G( x , y )

Gradient Descent Optimizer L too large ∆ G( x , y ) = f( x , y ) S = L ∙ G( x , y )

Gradient Descent Optimizer L too small ∆ G( x , y ) = f( x , y ) S = L ∙ G( x , y )

Gradient Descent Optimizer What’s wrong with this algorithm ?

Gradient Descent Optimizer S Units ? f(x, y) Units ? G(x, y) Units ? = millimeters = intensity / millimeters S = L ∙ G( x , y ) L Units ? = millimeters 2 / intensity

Gradient Descent Optimizer S = L ∙ G( x ) 1 1 f( x )

Gradient Descent Optimizer S = L ∙ G( x ) f( x ) S = Large in high gradients S = Small in low gradients 1 1

Gradient Descent Optimizer f( x ) S = L ∙ G( x ) Works great with this function Works badly with this function

Gradient Descent Variant Driving Safe !

Regular Step Gradient Descent f( x ) ^ S = D ∙ G( x ) If G changes direction Di = Di-1 / 2. 0 D 1 D 2 D 1

Regular Step Gradient Descent f( x ) ^ S = D ∙ G( x ) User Selects D 1 User Selects Dstop User Selects Gstop D 1 D 2 D 1

Optimizers are like a car Watch while you are driving !

Watch over your optimizer Example: Optimizer registering an image with itself starting at (-15 mm, -25 mm)

Mean Squared Differences Step Length = 1. 0 mm

Mean Squared Differences Step Length = 2. 0 mm

Mean Squared Differences Step Length = 5. 0 mm

Mean Squared Differences Step Length = 10. 0 mm

Mean Squared Differences Step Length = 20. 0 mm

Mean Squared Differences Step Length = 40. 0 mm

Watch over your optimizer Example: Optimizer registering an image shifted by ( -15 mm, -25 mm) The optimizer starts at (0 mm, 0 mm)

Mean Squared Differences Step Length = 1. 0 mm

Mean Squared Differences Step Length = 2. 0 mm

Mean Squared Differences Step Length = 5. 0 mm

Mean Squared Differences Step Length = 10. 0 mm

Mean Squared Differences Step Length = 20. 0 mm

Mean Squared Differences Step Length = 40. 0 mm

Multi - Modality

Multi-Modality Registration Fixed Image Moving Image Registered Moving Image

Multiple Image Modalities Number of pairs

Multiple Image Modalities More possible pairs

Mutual Information

Mutual Information

Exercise 21

Image Registration Fixed Image Moving Image Registered Moving Image

Image Registration #include ”itk. Image. Registration. Method. h” #include ”itk. Translation. Transform. h” #include ”itk. Mattes. Mutual. Information. Image. To. Image. Metric. h” #include ”itk. Linear. Interpolate. Image. Function. h” #include ”itk. Regular. Step. Gradient. Descent. Optimizer. h” #include ”itk. Image. File. Reader. h” #include ”itk. Image. File. Writer. h” #include ”itk. Resample. Image. Filter. h”

Image Registration const unsigned int Dimension = 2; typedef unsigned char Pixel. Type; typedef typedef itk: : Image< Pixel. Type , Dimension > itk: : Translation. Transform< double, Dimension > itk: : Regular. Step. Gradient. Descent. Optimizer itk: : Linear. Interpolate. Image. Function< Moving. Image. Type , double > Fixed. Image. Type; Moving. Image. Type; Transform. Type; Optimizer. Type; Interpolator. Type; typedef itk: : Mattes. Mutual. Information. Image. To. Image. Metric< Fixed. Image. Type , Moving. Image. Type > Metric. Type; typedef itk: : Image. Registration. Method< Fixed. Image. Type , Moving. Image. Type > Registration. Type;

Image Registration Transform. Type: : Pointer Optimizer. Type: : Pointer Interpolator. Type: : Pointer Metric. Type: : Pointer Registration. Type: : Pointer transform optimizer interpolator metric registrator->Set. Transform( registrator->Set. Optimizer( registrator->Set. Interpolator( registrator->Set. Metric( = Transform. Type: : New(); = Optimizer. Type: : New(); = Interpolator. Type: : New(); = Metric. Type: : New(); = Registration. Type: : New(); transform ); optimizer ); interpolator ); metric ); registrator->Set. Fixed. Image( fixed. Image. Reader->Get. Output() ); registrator->Set. Moving. Image( moving. Image. Reader->Get. Output() );

Image Registration registrator->Set. Fixed. Image. Region( fixed. Image. Reader->Get. Output()->Get. Buffered. Region() ); typedef Registration. Type: : Parameters. Type; transform->Set. Identity(); registrator->Set. Initial. Transform. Parameters( transform->Get. Parameters() ); optimizer->Set. Maximum. Step. Length( 4. 00 ); optimizer->Set. Minimum. Step. Length( 0. 05 ); optimizer->Set. Number. Of. Iterations( 200 ); optimizer->Maximize. On();

Image Registration metric->Set. Number. Of. Histogram. Bins( 20 ); metric->Set. Number. Of. Spatial. Samples( 10000 ); // Metric specific try { registrator->Start. Registration (); } catch( itk: : Exception. Object & excp ) { std: : cerr << “Error in registration” << std: : endl; std: : cerr << excp << std: : endl; } transform->Set. Parameters( registrator->Get. Last. Transform. Parameters() );

Image Registration typedef itk: : Resample. Image. Filter< Fixed. Image. Type , Moving. Image. Type > Resampler. Type; Resampler. Type : : Pointer resampler = Resampler. Type: : New(); resampler->Set. Transform ( transform ); resampler->Set. Input( moving. Image. Reader->Get. Output() ); Fixed. Image. Type : : Pointer fixed. Image = fixed. Image. Reader->Get. Output(); resampler->Set. Origin( fixed. Image->Get. Origin() ); resampler->Set. Spacing( fixed. Image->Get. Spacing() ); resampler->Set. Size( fixed. Image->Get. Largest. Possible. Region()->Get. Size() ); resampler->Update();

Rotation – Translation Parameter Scaling

Exercise 22

Image Registration Fixed Image Moving Image Registered Moving Image

Image Registration #include ”itk. Image. Registration. Method. h” #include ”itk. Centered. Rigid 2 DTransform. h” #include ”itk. Mean. Squares. Image. To. Image. Metric. h” #include ”itk. Linear. Interpolate. Image. Function. h” #include ”itk. Regular. Step. Gradient. Descent. Optimizer. h” #include ”itk. Centered. Transform. Initializer. h” #include ”itk. Image. File. Reader. h” #include ”itk. Image. File. Writer. h” #include ”itk. Resample. Image. Filter. h”

Image Registration const unsigned int Dimension = 2; typedef unsigned char Pixel. Type; typedef itk: : Image< Pixel. Type , Dimension > Fixed. Image. Type; Moving. Image. Type; typedef itk: : Centered. Rigid 2 DTransform< double > Transform. Type; typedef itk: : Centered. Transform. Initializer< Transform. Type , Fixed. Image. Type , Moving. Image. Type > Initializer. Type;

Image Registration Transform. Type: : Pointer Initializer. Type: : Pointer Optimizer. Type: : Pointer Interpolator. Type: : Pointer Metric. Type: : Pointer Registration. Type: : Pointer transform initializer optimizer interpolator metric registrator->Set. Transform( registrator->Set. Optimizer( registrator->Set. Interpolator( registrator->Set. Metric( = Transform. Type: : New(); = Initializer. Type: : New(); = Optimizer. Type: : New(); = Interpolator. Type: : New(); = Metric. Type: : New(); = Registration. Type: : New(); transform ); optimizer ); interpolator ); metric ); registrator->Set. Fixed. Image( fixed. Image. Reader->Get. Output() ); registrator->Set. Moving. Image( moving. Image. Reader->Get. Output() );

Image Registration registrator->Set. Fixed. Image. Region( fixed. Image. Reader->Get. Output()->Get. Buffered. Region() ); initializer->Set. Transform ( transform ); initializer->Set. Fixed. Image( fixed. Image. Reader->Get. Output() ); initializer->Set. Moving. Image( moving. Image. Reader->Get. Output() ); initializer->Moments. On(); initializer->Initialize. Transform(); registrator->Set. Initial. Transform. Parameters( transform->Get. Parameters() );

Image Registration typedef Optimizer. Type: : Scale. Type Optimizer. Scales. Type; Optimizer. Scales. Type optimizer. Scales( optimizer->Set. Maximum. Step. Length() ); const double translation. Scale = 1. 0 / 1000. 0 ; optimizer. Scales[ 0 ] = 1. 0; optimizer. Scales[ 1 ] = translation. Scale; optimizer. Scales[ 2 ] = translation. Scale; optimizer. Scales[ 3 ] = translation. Scale; optimizer. Scales[ 4 ] = translation. Scale; optimizer->Set. Scales( optimizer. Scales );

Image Registration try { registrator->Start. Registration (); } catch( itk: : Exception. Object & excp ) { std: : cerr << “Error in registration” << std: : endl; std: : cerr << excp << std: : endl; } transform->Set. Parameters( registrator->Get. Last. Transform. Parameters() );

Image Registration typedef itk: : Resample. Image. Filter< Fixed. Image. Type , Moving. Image. Type > Resampler. Type; Resampler. Type : : Pointer resampler = Resampler. Type: : New(); resampler->Set. Transform ( transform ); resampler->Set. Input( moving. Image. Reader->Get. Output() ); Fixed. Image. Type : : Pointer fixed. Image = fixed. Image. Reader->Get. Output(); resampler->Set. Origin( fixed. Image->Get. Origin() ); resampler->Set. Spacing( fixed. Image->Get. Spacing() ); resampler->Set. Size( fixed. Image->Get. Largest. Possible. Region()->Get. Size() ); resampler->Update();

Multi - Resolution

Multi-Resolution Registration Framework Improve speed Accuracy Robustness

Multi-Resolution Registration Framework Registration Fixed Image Moving Image

Multi-Resolution Registration Framework Flexible framework allows change of Parameters Components between resolution levels

Enjoy ITK !