[Insight-users] Implementation Finite Element Method (FEM) in ITK

[brian avants]

jeremy

> Where is this actually happening in the code. Is there any kind of
> documentation about this.

these are determined by each element --- we use standard approaches.
itkFEMElementBase.h provides a generic implementation.

>
> When and how is the global linear system created?

once you set up a mesh and all of the associated parameters, then you
know the # of degrees of freedom.  this determines the size of the
matrix.  classic FEM tells you exactly how to map the elements into
the global matrix given their geometry.   so, your question is usually
answered by the first chapter or so in a FEM book.

within the FEMCrankNicolsonSolver, there is a function called
AssembleKandM that does this in a generic way, assuming both a K and M
matrix.

b

>
> Again, I really appreciate your input.
>
> -J
>
>
> On Wed, Feb 2, 2011 at 10:19 PM, brian avants <stnava at gmail.com> wrote:
>>
>> jeremy
>>
>> sorry, i was really responding to a question cristina asked some time
>> ago and it took me a bit to remember the relevant pieces of math to
>> properly answer her .
>>
>> > Are you able to explain the bigger picture?
>>
>> in fem registration, you want to "deform" an image slowly over time
>> until the external image forces reach equilibrium with the internal
>> body forces (the regularization provided by the physical model).  the
>> parameters of the model control all of the above interactions.
>>
>> > I have the reference and the target image.
>> > If no grid/mesh is specified a uniform mesh is created.
>>
>> yes, i think that's correct.
>>
>> > I assume there are no initial forces. But what happens then?
>>
>> there is one fixed boundary condition to make the problem well-posed.
>> otherwise, no forces.
>>
>> > How are the forces and deformations for every node in every element
>> > calculated?
>>
>> we use standard gaussian integration points to integrate the
>> derivatives of an image similarity over each element and distribute
>> these derivatives to the nodes of the elements.
>>
>> > Is it based on pixel data (histogram matching) close to the nodes?
>>
>> there are a few options but you've got the basic idea --- you're
>> trying to match the local intensity.
>>
>> > Based on the local energy at the element level, how is the global energy
>> > calculated?
>>
>> not sure exactly what you mean here, but the elements are placed into
>> the global linear system that models the whole problem in the physical
>> space.
>>
>> b.a.
>>
>> > These are the kind of questions I have, it be really awesome if you can
>> > help.
>> > Thanks,
>> > -J
>> >
>> >
>> > On Tue, Feb 1, 2011 at 8:17 PM, brian avants <stnava at gmail.com> wrote:
>> >>
>> >> hi jeremy , cristina
>> >>
>> >> here is a brief latex sketch of the method used in the
>> >> FEMCrankNicolsonSolver --
>> >>
>> >> \text{matrix form parabolic PDE, heat equation type} \\
>> >> M\frac{ \partial U}{\partial t} +  K U  = f \\
>> >> \frac{ M U_t - M  U_{t-1}}{\delta_t} +  K U  = f \\
>> >> \frac{ M U_t - M  U_{t-1} +  \delta_t K U }{\delta_t} = f \\
>> >> \text{ here i make the decision to represent $U$ and $f$ at $t$ rather
>> >> than $t-1$ i believe this is backward Euler}\\
>> >> ( M  + \delta_t K ) U_t  =  M  U_{t-1} + \delta_t f_t \\
>> >> \\
>> >> \text{other choices are possible, e.g. } \\
>> >> \text{ Crank-Nicolson discretization is based on averaging forward and
>> >> backward Euler } \\
>> >> ( M + \alpha \delta_t K ) U_t = (  M - ( 1- \alpha ) \delta_t  K )
>> >> U_{t-1} + f \\
>> >>  f  =  \delta_t ( \alpha f_t + ( 1 - \alpha ) f_{t-1} ) \\
>> >> \\
>> >> \text{if}~ \alpha = 1 \text{ backward Euler }\\
>> >> ( M +  \delta_t K ) U_t = M  U_{t-1} + \delta_t f \\
>> >> \\
>> >> \text{if}~ \alpha = 0 \text{ forward Euler }\\
>> >>  M U_t =(  M - \delta_t  K ) U_{t-1} + \delta_t f_{t-1}
>> >>
>> >> hopefully this is helpful.
>> >>
>> >> brian
>> >>
>> >>
>> >>
>> >>
>> >> On Tue, Feb 1, 2011 at 4:01 PM, Jeremy Bournesel
>> >> <jeremy.bournesel at gmail.com> wrote:
>> >> > Hi,
>> >> > I have a decent background regarding FEM in the engineering domain
>> >> > and
>> >> > have
>> >> > started looking at it for registering medical images as well.
>> >> > I ran one of the ITK deformable registration demos, read the Software
>> >> > Guide
>> >> > and the ITK Powerpoint Presentations, however there is still some
>> >> > magic
>> >> > going on that I don't completely get.
>> >> > From what I understand ITK is using a uniform grid/mesh (if a custom
>> >> > one
>> >> > wasn't supplied) and then calculates the deformation for every
>> >> > element
>> >> > at
>> >> > the nodal points (iteratively).
>> >> > I couldn't find any information on how the exact process is working
>> >> > (besides
>> >> > some high level slides of the type K U = F).
>> >> > It'd be awesome if someone can point me to the right document or
>> >> > explain
>> >> > it
>> >> > to me.
>> >> > Thanks,
>> >> > Jeremy
>> >> >
>> >> >
>> >> >
>> >> >
>> >> >
>> >> > _____________________________________
>> >> > Powered by www.kitware.com
>> >> >
>> >> > Visit other Kitware open-source projects at
>> >> > http://www.kitware.com/opensource/opensource.html
>> >> >
>> >> > Kitware offers ITK Training Courses, for more information visit:
>> >> > http://www.kitware.com/products/protraining.html
>> >> >
>> >> > Please keep messages on-topic and check the ITK FAQ at:
>> >> > http://www.itk.org/Wiki/ITK_FAQ
>> >> >
>> >> > Follow this link to subscribe/unsubscribe:
>> >> > http://www.itk.org/mailman/listinfo/insight-users
>> >> >
>> >> >
>> >
>> >
>
>