This is the first example


Get default toolbox configurtion

myElxConf = elxDefaultConfiguration();

Feel free to adjust the settings in

In those files, you may adjust:

You probably do not need to adjust the other fields which are filenames skeletons.

Test elastix version

[Status, ErrorMessage] = elxTestElastixVersion(myElxConf);
if ~Status,

Read the street color image and convert it to a gray image

I = imread('street1.jpg');
I = 0.2989*I(:,:,1) + 0.5870*I(:,:,2) + 0.1140*I(:,:,3);
I = single(I);
ISize = size(I);
MaxI = max(I(:));

Define the fixed image

Elastix expects images to be a structure StrDatax structure with two fields Data for the image value, and x for the points coordinates.

FixedImage.Data = I;
FixedImage.x{1} = 5*(0:ISize(1)-1);
FixedImage.x{2} = 5*(0:ISize(2)-1);
MeanFixedImage = mean(FixedImage.Data(:));

Define the moving image

The moving image is the fixed image translated by 5*(-35, -15)=(-175,-75). We also add noise to the moving image.

MovingImage = FixedImage;
MovingImage.Data = MeanFixedImage*ones(ISize);
MovingImage.Data(1:end-35,1:end-15) = I(36:end, 16:end);
MovingImage.Data = MovingImage.Data + 20*randn(ISize);

Plot the fixed and moving images

figure(1); colormap gray
subplot(1, 2, 1);
imagesc(FixedImage.x{[2 1]}, FixedImage.Data, [0 MaxI]);
axis image; title('Fixed image');
xlabel('x\{2\}'); ylabel('x\{1\}');
subplot(1, 2, 2);
imagesc(MovingImage.x{[2 1]}, MovingImage.Data, [0 MaxI]);
axis image; title('Moving image');
xlabel('x\{2\}'); ylabel('x\{1\}');

Choose the elastix parameters

We choose:

myElxParam = elxDefaultParameters('EulerTransform', 2);
myElxParam.Optimizer = 'QuasiNewtonLBFGS';
myElxParam.ImageSampler = 'Full';
NumberOfResolutions = 4;
myElxParam.Metric = 'AdvancedMeanSquares';
myElxParam.NumberOfResolutions = NumberOfResolutions;
myElxParam.BSplineInterpolationOrder = 1;
myElxParam.FixedImagePyramid = 'FixedRecursiveImagePyramid';
myElxParam.MovingImagePyramid = 'MovingRecursiveImagePyramid';
myElxParam.DefaultPixelValue = MeanFixedImage;

Remove useless default configuration fields

myElxParam = rmfield(myElxParam, {'SampleRegionSize',...
  'NewSamplesEveryIteration', ...

Run Elastix

[strRegMoving,Transforms,log,success,message] = elxElastix(myElxConf, ...
  myElxParam, FixedImage, MovingImage);

disp(['success : ' num2str(success)])
if ~success
  disp(sprintf('  elastix message: %s\n',message{:}))
success : 1

The rigid transform found by Elastix

ans =


Transforms is a cell because all the available transforms have not the same fields. And when you specified a composed transform, elastix return several transforms.

Plot the registerd moving image and the metric

figure(2); colormap gray
subplot(1, 2, 1)
imagesc(strRegMoving.x{[2 1]}, strRegMoving.Data); axis image;
title('Registered moving image');
xlabel('x\{2\}'); ylabel('x\{1\}');
subplot(1, 2, 2);
imagesc(strRegMoving.x{[2 1]}, strRegMoving.Data - FixedImage.Data)
axis image
title('Difference between the registered image and the fixed image');
xlabel('x\{2\}'); ylabel('x\{1\}');

With a quasi-Newton optimiser, the number of iterations is not fixed and varies from one resolution to another.

NbIter = ones([1 NumberOfResolutions]);
for cpt = 1:NumberOfResolutions
  NbIter(cpt) = numel(log.ParameterFile(1).Resolution(cpt).Metric);
Legend = cell(1, NumberOfResolutions);
Metric = NaN(max(NbIter), NumberOfResolutions);
for cpt = 1:NumberOfResolutions
  Legend{cpt} = sprintf('Resolution %i', cpt-1);
  Metric(1:NbIter(cpt), cpt) = log.ParameterFile(1).Resolution(cpt).Metric;

Do the real plotting

xlabel('Iteration #');

log structure

We explore the log fields.

log.Command: The system command ran by elxElastix. Useful for debugging.

ans =

LD_PRELOAD=/usr/lib/libstdc++.so.6 elastix -f /tmp/elxInput/fixed.000.mhd -m /tmp/elxInput/moving.000.mhd -out /tmp/elxOutput -p /tmp/elxInput/parameter.0.txt

ans =

     4   600

ans =


log.ParameterFile(1) is a rather complex structure built by parsing in the logfile the part devoted to the first elastix ParameterFile. Here we seek for one rigid transform so there is only one ParameterFile.

We find the list of parameters we specified.

ans = 

               FixedInternalImagePixelType: 'float'
              MovingInternalImagePixelType: 'float'
                       FixedImageDimension: 2
                      MovingImageDimension: 2
                                 Transform: 'EulerTransform'
                              Registration: [1x27 char]
                         FixedImagePyramid: [1x26 char]
                        MovingImagePyramid: [1x27 char]
                              ImageSampler: 'Full'
                              Interpolator: [1x19 char]
                                    Metric: [1x19 char]
                                 Optimizer: 'QuasiNewtonLBFGS'
                 BSplineInterpolationOrder: 1
                    HowToCombineTransforms: 'Compose'
                 MaximumNumberOfIterations: 300
                     UseRandomSampleRegion: 0
                      ResampleInterpolator: [1x24 char]
                                 Resampler: 'DefaultResampler'
                       NumberOfResolutions: 4
                                 ErodeMask: 1
            FinalBSplineInterpolationOrder: 1
                      ResultImagePixelType: 'float'
                         DefaultPixelValue: 101.6050
     WriteTransformParametersEachIteration: 0
    WriteTransformParametersEachResolution: 0
       WriteResultImageAfterEachResolution: 0
                       UseDirectionCosines: 1

Default parameters are those that elastix have chosen for you. In the logfile they appear as warning messages.

Take care Elastix may choese a parameter because the name of the field is incorrect. So it is worth exploring this parameter.

ans = 

          AutomaticTransformInitialization: 0
    AutomaticTransformInitializationMethod: 'GeometricalCenter'
                 AutomaticScalesEstimation: 0
              GenerateLineSearchIterations: 0

Unprocessed warnings: They are informative messages about your configuraiton.

ans =

WARNING: the fixed pyramid schedule is not fully specified!

ans =

  A default pyramid schedule is used.

ans =

WARNING: the moving pyramid schedule is not fully specified!

ans =

  A default pyramid schedule is used.

ans =


For each resolution we get info

ans = 

       DefaultParameters: [1x1 struct]
     UnprocessedWarnings: {}
               SrchDirNr: [7x1 double]
                LineItNr: [7x1 double]
                  Metric: [7x1 double]
              StepLength: [7x1 double]
            GradientNorm: [7x1 double]
           SearchDirNorm: [7x1 double]
             DirGradient: [7x1 double]
                   Phase: {7x1 cell}
                  Wolfe1: [7x1 logical]
                  Wolfe2: [7x1 logical]
    LinSrchStopCondition: {7x1 cell}
                  Timems: [7x1 double]

ans = 

               NewSamplesEveryIteration: 0
                   ShowExactMetricValue: 0
                   CheckNumberOfSamples: 1
                       UseNormalization: 0
          NumberOfSamplesForSelfHessian: 100000
              SelfHessianSmoothingSigma: 1
                  SelfHessianNoiseRange: 1
    MaximumNumberOfLineSearchIterations: 20
                             StepLength: 1
               LineSearchValueTolerance: 1.0000e-04
            LineSearchGradientTolerance: 0.9000
             GradientMagnitudeTolerance: '1e-06'
                    LBFGSUpdateAccuracy: 5
                StopIfWolfeNotSatisfied: 1

Apply the transformation to a set of points

We define four points from their indices. Here the convention is the one of elastix not Matlab so one point have indices (0, 0). We could also have defined them with their coordinates.

PointSet.Indices = [0 0; 1 1; 10 1; 1 10];
[Data, Log, Success, Message]=elxTransformix(myElxConf, Transforms, ...
  'Deformation', PointSet);
Data = 

          InputPointSet: [1x1 struct]
         OutputPointSet: [1x1 struct]
    DeformationPointSet: [4x2 double]

ans =

 -173.5641  -76.4852
 -168.5624  -71.4869
 -123.5624  -71.5020
 -168.5473  -26.4869


Copyright (C) CNRS and Riverside Research Contributors: Alain CORON, Jonathan MAMOU (2010)

alain.coron@upmc.fr, JMamou@riversideresearch.org

This software is a computer program whose purpose is to effectively register images within Matlab (http://www.mathworks.com) with elastix (http://elastix.isi.uu.nl/), an open-source image-registration software.

This software was supported in part by NIH Grant CA100183, the Riverside Research Biomedical Engineering Research Fund, and CNRS.

This software is governed by the CeCILL-B license under French law and abiding by the rules of distribution of free software. You can use, modify and/ or redistribute the software under the terms of the CeCILL-B license as circulated by CEA, CNRS and INRIA at the following URL "http://www.cecill.info".

As a counterpart to the access to the source code and rights to copy, modify and redistribute granted by the license, users are provided only with a limited warranty and the software's author, the holder of the economic rights, and the successive licensors have only limited liability.

In this respect, the user's attention is drawn to the risks associated with loading, using, modifying and/or developing or reproducing the software by the user in light of its specific status of free software, that may mean that it is complicated to manipulate, and that also therefore means that it is reserved for developers and experienced professionals having in-depth computer knowledge. Users are therefore encouraged to load and test the software's suitability as regards their requirements in conditions enabling the security of their systems and/or data to be ensured and, more generally, to use and operate it in the same conditions as regards security.

The fact that you are presently reading this means that you have had knowledge of the CeCILL-B license and that you accept its terms.

$Id: elxExampleStreet.m 3 2012-05-25 19:37:07Z coron $