[Insight-users] Deformable Registration
Luca Nicotra
nicotra at cli.di.unipi.it
Sun Jun 6 13:11:04 EDT 2004
Thank you Luis,
the double parameters were just an experiment (I forgot to change the
file back)...but you are right, I forgot the
X->SetConfigFileName(paramname).
Cheers
-luca-
On Sun, 2004-06-06 at 11:33, Luis Ibanez wrote:
> 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