[ITK] On the deformable model
Nicolás Barnafi
nabw91 at gmail.com
Fri Mar 24 17:19:41 EDT 2017
Hi Matt,
I think I haven't explained myself correctly. At the end of the
presentation there is a variational formulation of the problem, different
to the one shown at the Optical Flow section. In that one, the problem reads
D[R,T;u] + a S[u],
which you can derivate to get a nonlinear system of equations
Au = f(u).
This is an SSD metric with an elastic penalization (assuming S is such a
functional) and as far as I understand it should be the default metric as
stated in the software guide. This setting has no temporal evolution, and
can be solved by a fixed point iterative scheme as
A u[k+1] = f(u[k])
or through a Newton-Raphson iterative scheme by looking for solutions of
F(u) := Au - f(u) = 0,
where you take a derivative and then iterate taking F(u[k+1]) = 0. Up to
this point, there is no time in any sense, which is confusing from the
example. When I read the code from DeformableRegistration1.h (the example I
mentioned before) and follow all the way what is happening, first the
metric is loaded
https://github.com/InsightSoftwareConsortium/ITK/blob/aecfac233e8815cdd0121fd2351dd3fef80d2e1b/Modules/Registration/FEM/include/itkFEMRegistrationFilter.hxx#L536-L538
and then the solver is called to iterate
https://github.com/InsightSoftwareConsortium/ITK/blob/aecfac233e8815cdd0121fd2351dd3fef80d2e1b/Modules/Registration/FEM/include/itkFEMRegistrationFilter.hxx#L547-L549
.
Now, if I go to check the kind of solver used, it is the CrankNicolson
Solver
https://github.com/InsightSoftwareConsortium/ITK/blob/aecfac233e8815cdd0121fd2351dd3fef80d2e1b/Modules/Registration/FEM/include/itkFEMRegistrationFilter.h#L145
and checking the documentation for it says explicitly it is used for time
dependent problems (
https://itk.org/Doxygen/html/classitk_1_1fem_1_1SolverCrankNicolson.html),
which the mentioned formulation is not. Now, in the variational
Registration filter (
https://itk.org/Doxygen/html/classitk_1_1VariationalRegistrationFilter.html),
it shows a way of solving this problem by creating an artificial time which
I am not sure if it is the one being used. For this, take the fixed point
scheme I mentioned, multiply by a time step and then use that the previous
solution is the steady state solution of the new problem:
u[k+1] + t A u[k+1] = t f(u[k]) + u[k].
This is also mentioned in Modersitzki's book (
http://www.oxfordscholarship.com/view/10.1093/acprof:oso/9780198528418.001.0001/acprof-9780198528418-chapter-9)
as a way of stabilizing the scheme, but also mentions the Optical Flow
formulation (no details). It is very tempting to think that the optical
flow formulation is being used, because it is mentioned in the elastic
registration ppt (the one we have been referring to) and because if I go to
the metrics you gave me, the one that uses the sum of squared differences
is the first one, which states in the documentation (
https://itk.org/Doxygen/html/classitk_1_1MeanSquareRegistrationFunction.html)
that it is actually used for the demons algorithm, again a very time
dependent way of approaching the problem, as it heavily depends on the
Optical Flow equations. Also, the solver I mentioned defines a timestep in
between, which again enforces the time idea.
What I am trying to do is use FEniCS for an alternative approach which
can't be implemented in ITK's current version, and as ITK is very stable
and robust, I would like to use it to see what is the correct solution, or
if it is really a problem of my data (although I have been using the brain
example for a week now). Right now I use the fixed point scheme, but my
solutions aren't stable at all, and even if they started being stable using
the Optical Flow formulation, I wouldn't be quite sure if that is what is
actually being used by ITK. I hope this clears the picture, my problem
isn't about mixing the concepts of iterations and timesteps, but about the
initial mathematical model used by the default metric, which shouldn't have
time steps at all. Anyway your code references have been great starting
points to read other sections, so thanks for everything so far.
Best regards
On 24 March 2017 at 16:37, Matt McCormick <matt.mccormick at kitware.com>
wrote:
> Hi Nicolás,
>
> Later on in the presentation it mention the iterative approach; some
> details are here:
>
> https://github.com/InsightSoftwareConsortium/ITK/blob/
> aecfac233e8815cdd0121fd2351dd3fef80d2e1b/Modules/Registration/FEM/include/
> itkFEMRegistrationFilter.hxx#L222-L255
>
> Within a FEM system solution, the similarity metric is not changing,
> but the metric is updated during iterations.
>
> HTH,
> Matt
>
> On Fri, Mar 24, 2017 at 3:20 PM, Nicolás Barnafi <nabw91 at gmail.com> wrote:
> > Hi Matt,
> >
> > Yes, and in the example the they don't set anything, so the metric is by
> > default the SSD. The part I don't see is, SSD does not include any kind
> of
> > evolution in time. With this in mind, does the solver assume, when we set
> > the use of a mass_matrix, that it has to use optical flow instead of
> simply
> > assessing the error given by the warping? If this is true, it would be
> great
> > to see the original optimization problem, or at least the EL equations
> and
> > how they are penalized. I am very intrigued because my (very) simple
> python
> > implementation uses no time evolution and a plain fixed point scheme,
> and no
> > matter how much I beg, it does not converge to anything that lungs like a
> > registration. Thanks again.
> >
> > Best regards
> >
> > On 24 March 2017 at 15:55, Matt McCormick <matt.mccormick at kitware.com>
> > wrote:
> >>
> >> Hi Nicolás,
> >>
> >> On slide 5, E_D and E_S are referring to the optical flow registration.
> >>
> >> On slide 6 is an objective function for variational registration, and
> >> FEM registration uses the integral of the deformation (linear elastic
> >> internal strain energy) as a regularization term.
> >>
> >> The similarity term can vary:
> >>
> >>
> >> https://github.com/InsightSoftwareConsortium/ITK/blob/
> aecfac233e8815cdd0121fd2351dd3fef80d2e1b/Modules/Registration/FEM/include/
> itkFEMRegistrationFilter.hxx#L222-L255
> >>
> >> HTH,
> >> Matt
> >>
> >> On Fri, Mar 24, 2017 at 11:39 AM, Nicolás Barnafi <nabw91 at gmail.com>
> >> wrote:
> >> > Hi Matt,
> >> >
> >> > Thanks the quick answer and welcome :) . That is indeed the ppt. There
> >> > they
> >> > define two energies and then write an integral in terms of some
> >> > "similarity"
> >> > which wasn't really defined. Maybe it is the E_D term? Hopefully the
> >> > same
> >> > applies for smoothness (E_S), but I'm not sure. In the code I had
> >> > trouble
> >> > because RunRegistration uses MultiResSolve, and there they
> instantiate a
> >> > solver which uses Update(). I couldn't find Update in either .h or
> .hxx
> >> > files, and apparently update is redirected (somehow) to AssembleKandM.
> >> > The
> >> > Doxygen page for the CrankNicolson solver shows a formulation of the
> >> > method
> >> > but not for what problem, and that leaves other questions without
> >> > answer,
> >> > such as the boundary conditions, elements used for the FEM
> formulation,
> >> > etc.
> >> > Thanks so much for your time.
> >> >
> >> > Best regards
> >> >
> >> > On 23 March 2017 at 18:34, Matt McCormick <matt.mccormick at kitware.com
> >
> >> > wrote:
> >> >>
> >> >> Hi Nicolás,
> >> >>
> >> >> Welcome to ITK!
> >> >>
> >> >>
> >> >> Is this the PowerPoint?
> >> >>
> >> >> https://itk.org/CourseWare/Training/NonRigidRegistrationMethods
> >> >>
> >> >> Note that DeformableRegistration1.cxx (FEM registration) is discussed
> >> >> after the optical flow example.
> >> >>
> >> >>
> >> >> Some more information on the FEMRegistrationFilter can be found in
> its
> >> >> class documentation:
> >> >>
> >> >>
> >> >>
> >> >> https://itk.org/Doxygen/html/classitk_1_1fem_1_
> 1FEMRegistrationFilter.html
> >> >>
> >> >> Multiple metrics are available:
> >> >>
> >> >>
> >> >>
> >> >> https://itk.org/Doxygen/html/classitk_1_1fem_1_
> 1FEMRegistrationFilter.html#acc636f4752e592cd780503a5fbfeba82
> >> >>
> >> >>
> >> >> In general, "the code reveals all" and should be given the highest
> >> >> amount of trust.
> >> >>
> >> >>
> >> >> Hope this helps,
> >> >> Matt
> >> >>
> >> >> On Thu, Mar 23, 2017 at 4:32 PM, Nicolás Barnafi <nabw91 at gmail.com>
> >> >> wrote:
> >> >> > Hi everyone,
> >> >> >
> >> >> > I have been trying to understand exactly what is happening in the
> >> >> > DeformableRegistration1.h file without much success. The problem
> is:
> >> >> > according to the software guide, the problem being solved comes
> from
> >> >> > the
> >> >> > variational problem given by
> >> >> >
> >> >> > min D[image1, image2; u] + S[u]
> >> >> >
> >> >> > where D is just the SSD (or the L2 norm of (Im1 - Im2 o phi), with
> >> >> > phi
> >> >> > the
> >> >> > unknown displacement field) and S is a linear elastic potential.
> From
> >> >> > here
> >> >> > you get the euler lagrange equations (asuming some unspecified
> >> >> > boundary
> >> >> > ocndition) and solve it using some semi implicit newton-raphson
> >> >> > scheme.
> >> >> > This
> >> >> > is where it starts getting blurry, because the ITK ppt on
> deformable
> >> >> > registration first shows an optical flow formulation, which would
> >> >> > mean
> >> >> > that
> >> >> > the SSD metric isn't really what is being used, and also if I dig
> >> >> > deeper
> >> >> > in
> >> >> > the code, I find actually a Crank-Nicolson scheme being used, which
> >> >> > really
> >> >> > implies some kind of temporality that really does not exist in the
> >> >> > variational formulation. The only hint I have found was in
> >> >> > Modersitzki's
> >> >> > book where a fixed point scheme is artificially stabilized:
> >> >> >
> >> >> > A(u[k+1]) = f_u[k]
> >> >> >
> >> >> > => u[k+1] + t A(u[k+1]) = t f_u[k] + u[k].
> >> >> >
> >> >> > I would want to know what is exactly happening in that example to
> be
> >> >> > able to
> >> >> > validate an example I implemented in python. Thanks for your time.
> >> >> >
> >> >> > Best regards
> >> >> >
> >> >> >
> >> >> >
> >> >> > --
> >> >> > Nicolás Alejandro Barnafi Wittwer
> >> >> >
> >> >> > _______________________________________________
> >> >> > Community mailing list
> >> >> > Community at itk.org
> >> >> > http://public.kitware.com/mailman/listinfo/community
> >> >> >
> >> >
> >> >
> >> >
> >> >
> >> > --
> >> > Nicolás Alejandro Barnafi Wittwer
> >
> >
> >
> >
> > --
> > Nicolás Alejandro Barnafi Wittwer
> _______________________________________________
> Community mailing list
> Community at itk.org
> http://public.kitware.com/mailman/listinfo/community
>
--
Nicolás Alejandro Barnafi Wittwer
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20170324/edbe7b91/attachment-0001.html>
More information about the Community
mailing list