[Insight-users] Deformable Registration
Luca Nicotra
nicotra at cli.di.unipi.it
Thu Jun 3 16:53:54 EDT 2004
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
More information about the Insight-users
mailing list