[Insight-users] Trouble with example deformableRegistration1.cxx

Markus Weigert m.weigert at fz-juelich.de
Tue Feb 7 07:09:51 EST 2006


Hi @all,

when I try to run the deformableRegistration1 example, it
crashes with the following message:

Microsoft Windows XP [Version 5.1.2600]
(C) Copyright 1985-2001 Microsoft Corp.

C:\DFG-Projekt\RIP\build\RelWithDebInfo>registrationImagingPlattform 
params.txt
Reading config file...params.txt
Example configured. E 0 rho 0
 reading moving
 reading fixed
 beginning level 0
 scaling 1
 scaling 1
 scaling 1
 ElementsPerDim 512 512 52
File  not found!
 applying loads

no landmark file specified.
 num of LM loads 0
 landmarks done
 allocating deformation field
 load sizes [512, 512, 52]  image [512, 512, 52]
 load sizes [512, 512, 52]  image [512, 512, 52]
 energy 0

This application has requested the Runtime to terminate it in an unusual 
way.
Please contact the application's support team for more information.


I use the following configurations file:

% Note: the paths in the parameters assume you have the traditional
% ITK file hierarchy as shown below:
%
% ITK/Insight/Examples/Registration/DeformableRegistration1.cxx
% ITK/Insight/Examples/Data/RatLungSlice*
% ITK/Insight-Bin/bin/DeformableRegistration1
%
% ---------------------------------------------------------
% Parameters for the single- or multi-resolution techniques
% ---------------------------------------------------------
1 % Number of levels in the multi-res pyramid (1 = single-res)
1 % Highest level to use in the pyramid
1 1 % Scaling at lowest level of pyramid
8 % Number of pixels per element
1.e4 % Elasticity (E)
1.e4 % Density x capacity (RhoC)
1 % Image energy scaling (gamma) - sets gradient step size
2 % NumberOfIntegrationPoints
1 % WidthOfMetricRegion
20 % MaximumIterations
% -------------------------------
% Parameters for the registration
% -------------------------------
0 0.99 % Similarity metric (0=mean sq, 1 = ncc, 2=pattern int, 3=MI, 
5=demons)
1.0 % Alpha
0 % DescentDirection (1 = max, 0 = min)
0 % DoLineSearch (0=never, 1=always, 2=if needed)
1.e1 % TimeStep
0.5 % Landmark variance
0 % Employ regridding / enforce diffeomorphism ( >= 1 -> true)
% ----------------------------------
% Information about the image inputs
% ----------------------------------
512 % Nx (image x dimension)
512 % Ny (image y dimension)
52  % Nz (image z dimension - not used if 2D)
rigid_reg_CT1.hdr % ReferenceFileName
CT.hdr % TargetFileName

% -------------------------------------------------------------------
% The actions below depend on the values of the flags preceding them.
% For example, to write out the displacement fields, you have to set
% the value of WriteDisplacementField to 1.
% -------------------------------------------------------------------
0 % UseLandmarks? - read the file name below if this is true
- % LandmarkFileName
./RatLung_result % ResultsFileName (prefix only)
0 % WriteDisplacementField?
./RatLung_disp % DisplacementsFileName (prefix only)
0 % ReadMeshFile?
- % MeshFileName
END

I also changed the code of the example a little bit to work in 3D (see 
the following part of
the message).
I hope that I didn't do a mistake here.
Help would be very appreciated to find the cause of this error.
I use ITK Version 2.4.0 and compile in release with deb.inf. mode.

Best regards,
Markus
 




/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: DeformableRegistration1.cxx,v $
  Language:  C++
  Date:      $Date: 2005/11/19 16:31:50 $
  Version:   $Revision: 1.32 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif


#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itkRescaleIntensityImageFilter.h"
#include "itkHistogramMatchingImageFilter.h"

//  Software Guide : BeginLatex
//
// The finite element (FEM) library within the Insight Toolkit can be
// used to solve deformable image registration problems.  The first step in
// implementing a FEM-based registration is to include the appropriate
// header files.
//
//  \index{Registration!Finite Element-Based}
//
//  Software Guide : EndLatex


// Software Guide : BeginCodeSnippet
#include "itkFEM.h"
#include "itkFEMRegistrationFilter.h"
      
// Software Guide : EndCodeSnippet

//#include "itkFEMFiniteDifferenceFunctionLoad.h"


//  Software Guide : BeginLatex
//
//  Next, we use \code{typedef}s to instantiate all necessary classes.  We
//  define the image and element types we plan to use to solve a
//  two-dimensional registration problem.  We define multiple element
//  types so that they can be used without recompiling the code.
//
//  Software Guide : EndLatex


//  Software Guide : BeginCodeSnippet
typedef itk::Image<unsigned char, 2>                       fileImageType;
typedef itk::Image<float, 2>                               ImageType;
typedef itk::fem::Element2DC0LinearQuadrilateralMembrane   ElementType;
typedef itk::fem::Element2DC0LinearTriangularMembrane      ElementType2;
//  Software Guide : EndCodeSnippet


//  Software Guide : BeginLatex
//
//  Note that in order to solve a three-dimensional registration
//  problem, we would simply define 3D image and element types in lieu
//  of those above.  The following declarations could be used for a 3D
//  problem:
//
//  SoftwareGuide : EndLatex


//  SoftwareGuide : BeginCodeSnippet
typedef itk::Image<short, 3>                            fileImage3DType;
typedef itk::Image<short, 3>                            Image3DType;
typedef itk::fem::Element3DC0LinearHexahedronMembrane   Element3DType;
typedef itk::fem::Element3DC0LinearTetrahedronMembrane  Element3DType2;
//  Software Guide : EndCodeSnippet


//  Software Guide : BeginLatex
// 
//  Here, we instantiate the load types and explicitly template the
//  load implementation type.  We also define visitors that allow the
//  elements and loads to communicate with one another. 
//
//  Software Guide : EndLatex


//typedef itk::fem::ImageMetricLoad<ImageType,ImageType>     ImageLoadType;

//  Software Guide : BeginCodeSnippet

typedef itk::fem::FiniteDifferenceFunctionLoad<Image3DType,Image3DType> 
ImageLoadType;
template class itk::fem::ImageMetricLoadImplementation<ImageLoadType>;

typedef Element3DType::LoadImplementationFunctionPointer     LoadImpFP;
typedef Element3DType::LoadType                              
ElementLoadType;

typedef Element3DType2::LoadImplementationFunctionPointer    LoadImpFP2;
typedef Element3DType2::LoadType                             
ElementLoadType2;

typedef itk::fem::VisitorDispatcher<Element3DType,ElementLoadType, 
LoadImpFP>  
                                                           DispatcherType;

typedef itk::fem::VisitorDispatcher<Element3DType2,ElementLoadType2, 
LoadImpFP2>  
                                                           DispatcherType2;
//  Software Guide : EndCodeSnippet


//  Software Guide : BeginLatex
//
//  Once all the necessary components have been instantiated, we can
//  instantiate the \doxygen{FEMRegistrationFilter}, which depends on the
//  image input and output types.
//
//  Software Guide : EndLatex


//  Software Guide : BeginCodeSnippet
typedef itk::fem::FEMRegistrationFilter<Image3DType,Image3DType> 
RegistrationType;
//  Software Guide : EndCodeSnippet


int main(int argc, char *argv[])
{
  char *paramname;
  if ( argc < 2 )
    {
    std::cout << "Parameter file name missing" << std::endl;
    std::cout << "Usage: " << argv[0] << " param.file" << std::endl;
    return -1;
    }
  else
    {
    paramname=argv[1];
    }


//  Software Guide : BeginLatex
// 
//  The \doxygen{fem::ImageMetricLoad} must be registered before it
//  can be used correctly with a particular element type.  An example
//  of this is shown below for ElementType.  Similar
//  definitions are required for all other defined element types.
//
//  Software Guide : EndLatex
 
  // Register the correct load implementation with the element-typed 
visitor dispatcher.
  {
//  Software Guide : BeginCodeSnippet
  Element3DType::LoadImplementationFunctionPointer fp =
    
&itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetricLoad;
  DispatcherType::RegisterVisitor((ImageLoadType*)0,fp);
//  Software Guide : EndCodeSnippet 
  }
  {
  Element3DType2::LoadImplementationFunctionPointer fp =
    
&itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetricLoad;
  DispatcherType2::RegisterVisitor((ImageLoadType*)0,fp);
  }


//  Software Guide : BeginLatex
// 
//  In order to begin the registration, we declare an instance of the
//  FEMRegistrationFilter.  For simplicity, we will call
//  it \code{registrationFilter}.
//
//  Software Guide : EndLatex

//  Software Guide : BeginCodeSnippet
  RegistrationType::Pointer registrationFilter = RegistrationType::New();
//  Software Guide : EndCodeSnippet


//  Software Guide : BeginLatex
//
//  Next, we call \code{registrationFilter->SetConfigFileName()} to read 
the parameter
//  file containing information we need to set up the registration
//  filter (image files, image sizes, etc.).  A sample parameter file is 
shown at the end of this
//  section, and the individual components are labeled. 
//
//  Software Guide : EndLatex


  // Attempt to read the parameter file, and exit if an error occurs
  registrationFilter->SetConfigFileName(paramname);
  if ( !registrationFilter->ReadConfigFile(
           (registrationFilter->GetConfigFileName()).c_str() ) )
    {
    return -1;
    }
 
  // Read the image files
  typedef itk::ImageFileReader< fileImage3DType >      FileSourceType;
  typedef fileImageType::PixelType PixType;

  FileSourceType::Pointer movingfilter = FileSourceType::New();
  movingfilter->SetFileName( 
"rigid_reg_ct1.hdr"/*(registrationFilter->GetMovingFile()).c_str()*/ );
  FileSourceType::Pointer fixedfilter = FileSourceType::New();
  fixedfilter->SetFileName( 
"ct.hdr"/*(registrationFilter->GetFixedFile()).c_str()*/ );
  std::cout << " reading moving " << registrationFilter->GetMovingFile() 
<< std::endl;
  std::cout << " reading fixed " << registrationFilter->GetFixedFile() 
<< std::endl;
 

  try
    {
    movingfilter->Update();
    }
  catch( itk::ExceptionObject & e )
    {
    std::cerr << "Exception caught during reference file reading " << 
std::endl;
    std::cerr << e << std::endl;
    return -1;
    }
  try
    {
    fixedfilter->Update();
    }
  catch( itk::ExceptionObject & e )
    {
    std::cerr << "Exception caught during target file reading " << 
std::endl;
    std::cerr << e << std::endl;
    return -1;
    }
 

  // Rescale the image intensities so that they fall between 0 and 255
  typedef itk::RescaleIntensityImageFilter<fileImage3DType,Image3DType> 
FilterType;
  FilterType::Pointer movingrescalefilter = FilterType::New();
  FilterType::Pointer fixedrescalefilter = FilterType::New();

  movingrescalefilter->SetInput(movingfilter->GetOutput());
  fixedrescalefilter->SetInput(fixedfilter->GetOutput());

  const double desiredMinimum =  0.0;
  const double desiredMaximum =  255.0;

  movingrescalefilter->SetOutputMinimum( desiredMinimum );
  movingrescalefilter->SetOutputMaximum( desiredMaximum );
  movingrescalefilter->UpdateLargestPossibleRegion();
  fixedrescalefilter->SetOutputMinimum( desiredMinimum );
  fixedrescalefilter->SetOutputMaximum( desiredMaximum );
  fixedrescalefilter->UpdateLargestPossibleRegion();
 

  // Histogram match the images
  typedef itk::HistogramMatchingImageFilter<Image3DType,Image3DType> 
HEFilterType;
  HEFilterType::Pointer IntensityEqualizeFilter = HEFilterType::New();

  IntensityEqualizeFilter->SetReferenceImage( 
fixedrescalefilter->GetOutput() );
  IntensityEqualizeFilter->SetInput( movingrescalefilter->GetOutput() );
  IntensityEqualizeFilter->SetNumberOfHistogramLevels( 100);
  IntensityEqualizeFilter->SetNumberOfMatchPoints( 15);
  IntensityEqualizeFilter->ThresholdAtMeanIntensityOn();
  IntensityEqualizeFilter->Update();

  registrationFilter->SetFixedImage(fixedrescalefilter->GetOutput());
  registrationFilter->SetMovingImage(IntensityEqualizeFilter->GetOutput());


  itk::ImageFileWriter<Image3DType>::Pointer writer;
  writer = itk::ImageFileWriter<Image3DType>::New();
  std::string ofn="fixed.hdr";
  writer->SetFileName(ofn.c_str());
  writer->SetInput(registrationFilter->GetFixedImage() );
  writer->Write();

  ofn="moving.hdr";
  itk::ImageFileWriter<Image3DType>::Pointer writer2;
  writer2 =  itk::ImageFileWriter<Image3DType>::New();
  writer2->SetFileName(ofn.c_str());
  writer2->SetInput(registrationFilter->GetMovingImage() );
  writer2->Write();
 



//  Software Guide : BeginLatex
//
//  In order to initialize the mesh of elements, we must first create
//  ``dummy'' material and element objects and assign them to the
//  registration filter.  These objects are subsequently used to
//  either read a predefined mesh from a file or generate a mesh using
//  the software.  The values assigned to the fields within the
//  material object are arbitrary since they will be replaced with
//  those specified in the parameter file.  Similarly, the element
//  object will be replaced with those from the desired mesh.
//
//  Software Guide : EndLatex
 
//  Software Guide : BeginCodeSnippet
  // Create the material properties
  itk::fem::MaterialLinearElasticity::Pointer m;
  m = itk::fem::MaterialLinearElasticity::New();
  m->GN = 0;                  // Global number of the material
  m->E = registrationFilter->GetElasticity();  // Young's modulus -- 
used in the membrane
  m->A = 1.0;                 // Cross-sectional area
  m->h = 1.0;                 // Thickness
  m->I = 1.0;                 // Moment of inertia
  m->nu = 0.;                 // Poisson's ratio -- DONT CHOOSE 1.0!!
  m->RhoC = 1.0;              // Density
 
  // Create the element type
  Element3DType::Pointer e1=Element3DType::New();
  e1->m_mat=dynamic_cast<itk::fem::MaterialLinearElasticity*>( m );
  registrationFilter->SetElement(e1);
  registrationFilter->SetMaterial(m);
//  Software Guide : EndCodeSnippet


//  Software Guide : BeginLatex
//
//  Now we are ready to run the registration:
//
//  Software Guide : EndLatex

//  Software Guide : BeginCodeSnippet
  registrationFilter->RunRegistration();
//  Software Guide : EndCodeSnippet


//  Software Guide : BeginLatex
//
//  To output the image resulting from the registration, we can call
//  \code{WriteWarpedImage()}.  The image is written in floating point
//  format.
//
//  Software Guide : EndLatex

//  Software Guide : BeginCodeSnippet
  registrationFilter->WriteWarpedImage(
        (registrationFilter->GetResultsFileName()).c_str());
//  Software Guide : EndCodeSnippet

//  Software Guide : BeginLatex
//
//  We can also output the displacement fields resulting from the
//  registration, we can call \code{WriteDisplacementField()} with the
//  desired vector component as an argument.  For a $2D$ registration,
//  you would want to write out both the $x$ and $y$ displacements, and
//  this requires two calls to the aforementioned function.
//
//  Software Guide : EndLatex

//  Software Guide : BeginCodeSnippet
  if (registrationFilter->GetWriteDisplacements())
    {
    //registrationFilter->WriteDisplacementField(0);
    //registrationFilter->WriteDisplacementField(1);
    // If this were a 3D example, you might also want to call this line:
    // registrationFilter->WriteDisplacementField(2);

    // We can also write it as a multicomponent vector field
    registrationFilter->WriteDisplacementFieldMultiComponent();
    }
//  Software Guide : EndCodeSnippet

  //  This is a documented sample parameter file that can be used with
  //  this deformable registration example.
  //
  //  ../Data/FiniteElementRegistrationParameters1.txt
  //

  return 0;
}



-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20060207/ebc58efa/attachment-0001.html


More information about the Insight-users mailing list