[Insight-users] FEM Initialization

[Michael Hardisty]

Hello ITK users,
                          I am currently trying to deformambly register 
two 3D CT images to each other.  They are both of the same bone, one is 
just under load and the other is not.  I am registering in an attempt to 
look at the deformation field.  I am currently attempting to use the FEM 
frame work with ITK.  I have two questions in regards to usage of FEM.

1)  I am doing a 3D registration of two uCT images and am getting the 
following error.  Does anyone know what is causing this error?:
"ExceptionObject caught !
Requested region is (at least partially) outside the largest possible 
region."

- I think that this exception may have something to do with the spacing 
of  my images.  which have spacing: "0.0180441 X 0.0194318 X 
0.0192933".  Is it a requirement that the FEMRegistrationFilter have 
spacing 1 X 1 X 1.  I have tracked this error to the following code in 
"itkFEMRegistrationFilter.txx"

this->ComputeJacobian(1.,m_Field,2.5);

2)  I wanted to initialize the filter from within the code and not use 
the file IO that is showed in "DeformableRegistration1.cxx" because I am 
putting this routine into a viewer as a plugin.  Is there an example or 
documentation about how to properly initialize everything that is 
necessary to use the FEMRegistrationFilter.  Many of the member 
functions seem to provide the capability of initialization of 
parameters.  However some parameters specified in the parameter file do 
not seem to have equivalent member functions that are capable of 
defining the same parameters. 
It is not clear how to set:  image dimensions, or the 
lowestLevelScaling, which are set in the param file.  I have put the 
relevant bits of my code below:

thank you

-- 
Michael Hardisty
M.A.Sc Student
University of Toronto
Orthopaedic Biomechanics Laboratory
Sunnybrook & Women's College Health Sciences Centre



//ITK PipeLine Beginning

       

        // Register the correct load implementation with the 
element-typed visitor dispatcher.
        {
            //  Software Guide : BeginCodeSnippet
            ElementType::LoadImplementationFunctionPointer fp =
                
&itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetricLoad;
                DispatcherType::RegisterVisitor((ImageLoadType*)0,fp);

        }
        {
            ElementType2::LoadImplementationFunctionPointer fp2 =
                
&itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetricLoad;
                DispatcherType2::RegisterVisitor((ImageLoadType*)0,fp2);
        }
       


        RegistrationType::Pointer registrationFilter = 
RegistrationType::New();


        //Configuring the FEMRegistrationFilter using function calls to 
Set Methods etc.
 
        int numMultiResLevels =            1;
        int highestLevel =                1;
        int lowestLevelScaling[3] =        {1,1,1};
        int pixelsPerElement =            4;
        float elasticity =                10000;
        float densityCapacity =            10000;
        float density =                    4;
        float EnergyScale =                1; //Sets Gradient Step Size
        int numIntegrationPoints =        2;
        int widthMetricRegion =            1;
        int maxIterations =                20;
        int similarityMetric =            0; //(0=mean sq, 1 = ncc, 
2=pattern int, 3=MI, 5=demons)
        float alpha =                    1;
        int descentDirection =            0; //(1 = max, 0 = min)
        int doLineSearch =                0; //(0=never, 1=always, 2=if 
needed)
        float timeStep =                10;
        float landMarkVariance =        0.5;
        int reGridding =                0; //Employ regridding / enforce 
diffeomorphism ( >= 1 -> true)
        long imageDims[ImageDimension];
        int useLandmarks =                0;
        int readMesh =                    0;
        int writeDisplacementField =    0;
        float gamma =                    0.75;

        registrationFilter->SetNumLevels(numMultiResLevels);
        registrationFilter->SetMaxLevel(highestLevel);
        //don't know How to use lowestLevelScaling
        int level = 0;
        
registrationFilter->SetMeshPixelsPerElementAtEachResolution(pixelsPerElement,0);
        
registrationFilter->SetMeshPixelsPerElementAtEachResolution(pixelsPerElement,1);
        
registrationFilter->SetMeshPixelsPerElementAtEachResolution(pixelsPerElement,2);
        registrationFilter->SetElasticity(elasticity);
          registrationFilter->SetRho(density,level);
        registrationFilter->SetEnergyReductionFactor(EnergyScale);
        
registrationFilter->SetNumberOfIntegrationPoints(numIntegrationPoints);
        registrationFilter->SetWidthOfMetricRegion(widthMetricRegion,level);
        registrationFilter->SetMaximumIterations(maxIterations,level);
        registrationFilter->ChooseMetric(similarityMetric);
        registrationFilter->SetAlpha(alpha);


       
        if (descentDirection == 0)
        {
            registrationFilter->SetDescentDirectionMinimize();
        }
        else
        {
            registrationFilter->SetDescentDirectionMaximize();
        }
        registrationFilter->DoLineSearch(doLineSearch);
        registrationFilter->SetTimeStep(timeStep);
        //don't know how to use landmark variance
        registrationFilter->EmployRegridding(reGridding);
        //don't know how to use imagedims
        registrationFilter->UseLandmarks(useLandmarks);
        //don't know how to use readMesh
        //not using writeDisplacementField
        registrationFilter->SetGamma(gamma);
        defImageLoaderType* movingImage = defImageLoader;
        fixImageLoaderType* fixedImage = fixImageLoader;

        registrationFilter->SetMovingImage(defImageLoader->GetOutput());
        registrationFilter->SetFixedImage(fixImageLoader->GetOutput());
 
 




        // Create the material properties
        itk::fem::MaterialLinearElasticity::Pointer m;
        m = itk::fem::MaterialLinearElasticity::New();
        m->GN = 0;                  // Global number of the material
        m->E = registrationFilter->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.29;                 // 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 );
        registrationFilter->SetElement(e1);
        registrationFilter->SetMaterial(m);


        InterpolatorType::Pointer   interpolator    = 
InterpolatorType::New();
        registrationFilter->SetInterpolator(interpolator);

        VectorFieldType::Pointer initVectorField = VectorFieldType::New();

        registrationFilter->SetDeformationField(initVectorField);

        try
        {
                       
            registrationFilter->RunRegistration();
                       
        }
        catch( itk::ExceptionObject & err )
        {
            theMsg->printf("ExceptionObject caught !");
            theMsg->printf(err.GetDescription());
            return;
                        //return;
        }

        RegistrationType::FixedImageType::Pointer warpedImage = 
registrationFilter->GetWarpedImage();
        RegistrationType::FieldType::Pointer deformationField = 
registrationFilter->GetDeformationField();