[Insight-users] Deformable Registration
Luis Ibanez
luis.ibanez at kitware.com
Sun Jun 6 11:15:07 EDT 2004
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
>
More information about the Insight-users
mailing list