[Insight-users] Deformable Registration

Luis Ibanez luis.ibanez at kitware.com
Sun Jun 6 11:33:08 EDT 2004


Hi Luca,

Also...
Your parameters file contain unnecessary entries:


> % ---------------------------------------------------------
> % Parameters for the single- or multi-resolution techniques
> % ---------------------------------------------------------
> 2	% Number of levels in the multi-res pyramid (1 = single-res)
> 2	% Highest level to use in the pyramid
>  16 16 16	% Scaling at lowest level of pyramid
>  4 4		% Number of pixels per element
>  1.e4 1.e4 	        % Elasticity (E)
>  1.e4 1.e4 	        % Density x capacity (RhoC)
>  1 1		% Image energy scaling (gamma) - sets gradient step size
>  4 1 		% NumberOfIntegrationPoints
>  4 4  		% WidthOfMetricRegion
>  10 10 		% MaximumIterations


You are providing two values for elasticity (!?),
and for RhoC, Energy scaling, integration points, width of metric,
and maximum iterations (!?)



For a 3D case, the entries should be:

% ---------------------------------------------------------
% Parameters for the single- or multi-resolution techniques
% ---------------------------------------------------------
2	% Number of levels in the multi-res pyramid (1 = single-res)
2	% Highest level to use in the pyramid
  16 16 16	% Scaling at lowest level of pyramid
  4 		% Number of pixels per element
  1.e4  	        % Elasticity (E)
  1.e4 	        % Density x capacity (RhoC)
  1 		% Image energy scaling (gamma) - sets gradient step size
  4  		% NumberOfIntegrationPoints
  4   		% WidthOfMetricRegion
  10  		% MaximumIterations



Those values are all scalars.




Regards,


    Luis


------------------------
Luis Ibanez wrote:

> 
> Hi Luca,
> 
> You missed to give the filename of the parameters
> file to the registration filter.
> 
> Please add a line like:
> 
>    X->SetConfigFileName(paramname);
> 
> before you invoke:
> 
>    X->ReadConfigFile(strParamFile.data());
> 
> 
> 
> This is indeed a redundant operation... and it will
> get deprecated at some point in the following months.
> This filter is scheduled for refactoring because it
> doesn't fully adhere to the style of the toolkit.
> In particular, it shouldn't be doing I/O by itself.
> 
> 
> 
> Please let us know if you find any other problem,
> 
> 
>   Thanks,
> 
> 
>     Luis
> 
> 
> ------------------
> Luca Nicotra wrote:
> 
>> Hi,
>> I am trying to use the FEMRegistrationFilter with 3D images, but when I
>> call RunRegistration() I get this exception:
>>
>> itk::FEMExceptionLinearSystem (0x84ca2f8)
>> Location: "LinearSystemWrapperItpack::InitializeSolution"
>> File:
>> /home/luca/ITK/InsightToolkit-1.6.0/Code/Numerics/FEM/itkFEMLinearSystemWr 
>> apperItpack.cxx
>> Line: 164
>> Description: Error in linear system: System order not set
>>
>> I am probably not initializing some value, or maybe my config file is
>> wrong, I don't know. Does anyone know where I am wrong?
>> Thank you in advance.
>>
>> -luca-
>>
>>
>>
>> These are my code and my config file
>>
>> #include "ImspRegistrationTools.h"
>>
>> #include "itkImageFileReader.h" #include "itkImageFileWriter.h" 
>> #include "itkRescaleIntensityImageFilter.h"
>> #include "itkHistogramMatchingImageFilter.h"
>> #include "itkFEM.h"
>> #include "itkFEMRegistrationFilter.h"
>>
>> namespace ImspLib {
>> typedef itk::Image<float, 3>                            ImageType;
>> typedef itk::fem::Element3DC0LinearHexahedronMembrane   ElementType;
>> typedef itk::fem::Element3DC0LinearTetrahedronMembrane  ElementType2;
>>
>> typedef itk::fem::FiniteDifferenceFunctionLoad<ImageType,ImageType>
>> ImageLoadType;
>> template class itk::fem::ImageMetricLoadImplementation<ImageLoadType>;
>>
>> typedef ElementType::LoadImplementationFunctionPointer     LoadImpFP;
>> typedef ElementType::LoadType                             
>> ElementLoadType;
>>
>> typedef ElementType2::LoadImplementationFunctionPointer    LoadImpFP2;
>> typedef ElementType2::LoadType                            
>> ElementLoadType2;
>>
>> typedef itk::fem::VisitorDispatcher<ElementType,ElementLoadType,
>> LoadImpFP>                                                             
>> DispatcherType;
>>
>> typedef itk::fem::VisitorDispatcher<ElementType2,ElementLoadType2,
>> LoadImpFP2>   
>>                                                           
>> DispatcherType2;
>>
>> typedef itk::fem::FEMRegistrationFilter<ImageType,ImageType>
>> RegistrationType;
>>
>> template <class TPixelType>
>> typename RegistrationTools<TPixelType>::ImagePointer3DType
>> ImspLib::RegistrationTools<TPixelType>::deformableRegistration(ImagePointer3DType 
>> pMovingImage,         ImagePointer3DType pFixedImage, std::string 
>> strParamFile)
>> {
>>
>>   {
>>   ElementType::LoadImplementationFunctionPointer fp =    
>> &itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetricLoad; 
>>
>>   DispatcherType::RegisterVisitor((ImageLoadType*)0,fp);   }
>>   {
>>   ElementType2::LoadImplementationFunctionPointer fp =
>>    
>> &itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetricLoad; 
>>
>>   DispatcherType2::RegisterVisitor((ImageLoadType*)0,fp);
>>   }
>>
>>   RegistrationType::Pointer X = RegistrationType::New();     
>> X->ReadConfigFile(strParamFile.data());
>>  
>>
>>   typedef itk::RescaleIntensityImageFilter<Image3DType,ImageType>
>> FilterType;
>>   typename FilterType::Pointer movingrescalefilter = FilterType::New();
>>   typename FilterType::Pointer fixedrescalefilter = FilterType::New();
>>
>>   movingrescalefilter->SetInput(pMovingImage);
>>   fixedrescalefilter->SetInput(pFixedImage);
>>
>>   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<ImageType,ImageType>
>> HEFilterType;
>>   HEFilterType::Pointer IntensityEqualizeFilter = HEFilterType::New();
>>
>>   IntensityEqualizeFilter->SetReferenceImage(
>> fixedrescalefilter->GetOutput() );
>>   IntensityEqualizeFilter->SetInput( movingrescalefilter->GetOutput() );
>>   IntensityEqualizeFilter->SetNumberOfHistogramLevels( 100);
>>   IntensityEqualizeFilter->SetNumberOfMatchPoints( 15);
>>   IntensityEqualizeFilter->ThresholdAtMeanIntensityOn();
>>   IntensityEqualizeFilter->Update();
>>
>>   X->SetFixedImage(fixedrescalefilter->GetOutput());
>>   X->SetMovingImage(IntensityEqualizeFilter->GetOutput());
>>
>>
>>   itk::ImageFileWriter<ImageType>::Pointer writer;
>>   writer = itk::ImageFileWriter<ImageType>::New();
>>   std::string ofn="fixed.hdr";
>>   writer->SetFileName(ofn.c_str());
>>   writer->SetInput(X->GetFixedImage() );   writer->Write();
>>
>>   ofn="moving.hdr";
>>   itk::ImageFileWriter<ImageType>::Pointer writer2;
>>   writer2 =  itk::ImageFileWriter<ImageType>::New();
>>   writer2->SetFileName(ofn.c_str());
>>   writer2->SetInput(X->GetMovingImage() );   writer2->Write();
>>  
>>
>>
>>
>>   itk::fem::MaterialLinearElasticity::Pointer m;
>>   m = itk::fem::MaterialLinearElasticity::New();
>>   m->GN = 0;                  // Global number of the material
>>   m->E = X->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   ElementType::Pointer 
>> e1=ElementType::New();
>>   e1->m_mat=dynamic_cast<itk::fem::MaterialLinearElasticity*>( m );
>>   X->SetElement(e1);
>>   X->SetMaterial(m);
>>   try
>>   {
>>   X->RunRegistration();
>>   }
>>     catch( itk::ExceptionObject & e )
>>     {
>>     std::cerr << "Exception caught during RunRegistration() " <<
>> std::endl;
>>     std::cerr << e << std::endl;
>>     }
>>   X->WriteWarpedImage("warped");
>>     X->WriteDisplacementFieldMultiComponent();
>>
>>   ImageType::Pointer bla = X->GetOutput();
>>   ImageType::Pointer cra = X->GetWarpedImage();
>>     typedef itk::RescaleIntensityImageFilter<ImageType, Image3DType>
>> ReturnFilterType;
>>   typename ReturnFilterType::Pointer  returnFilter=
>> ReturnFilterType::New();
>>     returnFilter->SetInput(bla);
>>   returnFilter->SetOutputMinimum(0);
>>   returnFilter->SetOutputMaximum(255);   return 
>> returnFilter->GetOutput();
>> }
>>
>> }
>>
>>
>>
>> //===================================CONFIG FILE========================
>>
>>
>>
>> % Configuration file #1 for DeformableRegistration1.cxx
>> % % This example demonstrates the setup of a basic registration % 
>> problem that does NOT use multi-resolution strategies.  As a % result, 
>> only one value for the parameters between % (# of pixels per element) 
>> and (maximum iterations) is necessary.
>> % If you were using multi-resolution, you would have to specify
>> % values for those parameters at each level of the pyramid.  %
>> % 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
>> % ---------------------------------------------------------
>> 2    % Number of levels in the multi-res pyramid (1 = single-res)
>> 2    % Highest level to use in the pyramid
>>  16 16 16    % Scaling at lowest level of pyramid
>>  4 4        % Number of pixels per element
>>  1.e4 1.e4             % Elasticity (E)
>>  1.e4 1.e4             % Density x capacity (RhoC)
>>  1 1        % Image energy scaling (gamma) - sets gradient step size
>>  4 1         % NumberOfIntegrationPoints
>>  4 4          % WidthOfMetricRegion
>>  10 10         % 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)
>> 20    % Nz (image z dimension - not used if 2D)
>> not_used.nihil     % ReferenceFileName not_used.nihil  % 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
>> ./bla_result                       % ResultsFileName (prefix only)
>> 1       % WriteDisplacementField?
>> ./bla_disp                         % DisplacementsFileName (prefix only)
>> 0       % ReadMeshFile?
>> -                                      % MeshFileName
>> END
>> _______________________________________________
>> Insight-users mailing list
>> Insight-users at itk.org
>> http://www.itk.org/mailman/listinfo/insight-users
>>
> 
> 
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
> 





More information about the Insight-users mailing list