[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