<div dir="ltr">Hi Matt, <div><br></div><div>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</div><div><br></div><div>D[R,T;u] + a S[u],</div><div><br></div><div>which you can derivate to get a nonlinear system of equations </div><div><br></div><div>Au = f(u). </div><div><br></div><div>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</div><div><br></div><div>A u[k+1] = f(u[k])</div><div><br></div><div>or through a Newton-Raphson iterative scheme by looking for solutions of</div><div><br></div><div>F(u) := Au - f(u) = 0, </div><div><br></div><div>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</div><div><br></div><div><a href="https://github.com/InsightSoftwareConsortium/ITK/blob/aecfac233e8815cdd0121fd2351dd3fef80d2e1b/Modules/Registration/FEM/include/itkFEMRegistrationFilter.hxx#L536-L538">https://github.com/InsightSoftwareConsortium/ITK/blob/aecfac233e8815cdd0121fd2351dd3fef80d2e1b/Modules/Registration/FEM/include/itkFEMRegistrationFilter.hxx#L536-L538</a><br></div><div><br></div><div>and then the solver is called to iterate</div><div><br></div><div><a href="https://github.com/InsightSoftwareConsortium/ITK/blob/aecfac233e8815cdd0121fd2351dd3fef80d2e1b/Modules/Registration/FEM/include/itkFEMRegistrationFilter.hxx#L547-L549">https://github.com/InsightSoftwareConsortium/ITK/blob/aecfac233e8815cdd0121fd2351dd3fef80d2e1b/Modules/Registration/FEM/include/itkFEMRegistrationFilter.hxx#L547-L549</a>.<br></div><div><br></div><div>Now, if I go to check the kind of solver used, it is the CrankNicolson Solver<br><br><a href="https://github.com/InsightSoftwareConsortium/ITK/blob/aecfac233e8815cdd0121fd2351dd3fef80d2e1b/Modules/Registration/FEM/include/itkFEMRegistrationFilter.h#L145">https://github.com/InsightSoftwareConsortium/ITK/blob/aecfac233e8815cdd0121fd2351dd3fef80d2e1b/Modules/Registration/FEM/include/itkFEMRegistrationFilter.h#L145</a><br></div><div><br></div><div>and checking the documentation for it says explicitly it is used for time dependent problems (<a href="https://itk.org/Doxygen/html/classitk_1_1fem_1_1SolverCrankNicolson.html">https://itk.org/Doxygen/html/classitk_1_1fem_1_1SolverCrankNicolson.html</a>), which the mentioned formulation is not. Now, in the variational Registration filter (<a href="https://itk.org/Doxygen/html/classitk_1_1VariationalRegistrationFilter.html">https://itk.org/Doxygen/html/classitk_1_1VariationalRegistrationFilter.html</a>), 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:</div><div><br></div><div>u[k+1] + t A u[k+1] = t f(u[k]) + u[k].</div><div><br></div><div>This is also mentioned in Modersitzki's book (<a href="http://www.oxfordscholarship.com/view/10.1093/acprof:oso/9780198528418.001.0001/acprof-9780198528418-chapter-9">http://www.oxfordscholarship.com/view/10.1093/acprof:oso/9780198528418.001.0001/acprof-9780198528418-chapter-9</a>) 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 (<a href="https://itk.org/Doxygen/html/classitk_1_1MeanSquareRegistrationFunction.html">https://itk.org/Doxygen/html/classitk_1_1MeanSquareRegistrationFunction.html</a>) 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. </div><div><br></div><div>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. </div><div><br></div><div>Best regards</div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On 24 March 2017 at 16:37, Matt McCormick <span dir="ltr"><<a href="mailto:matt.mccormick@kitware.com" target="_blank">matt.mccormick@kitware.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi Nicolás,<br>
<br>
Later on in the presentation it mention the iterative approach; some<br>
details are here:<br>
<br>
<a href="https://github.com/InsightSoftwareConsortium/ITK/blob/aecfac233e8815cdd0121fd2351dd3fef80d2e1b/Modules/Registration/FEM/include/itkFEMRegistrationFilter.hxx#L222-L255" rel="noreferrer" target="_blank">https://github.com/<wbr>InsightSoftwareConsortium/ITK/<wbr>blob/<wbr>aecfac233e8815cdd0121fd2351dd3<wbr>fef80d2e1b/Modules/<wbr>Registration/FEM/include/<wbr>itkFEMRegistrationFilter.hxx#<wbr>L222-L255</a><br>
<br>
Within a FEM system solution, the similarity metric is not changing,<br>
but the metric is updated during iterations.<br>
<br>
HTH,<br>
Matt<br>
<div class="HOEnZb"><div class="h5"><br>
On Fri, Mar 24, 2017 at 3:20 PM, Nicolás Barnafi <<a href="mailto:nabw91@gmail.com">nabw91@gmail.com</a>> wrote:<br>
> Hi Matt,<br>
><br>
> Yes, and in the example the they don't set anything, so the metric is by<br>
> default the SSD. The part I don't see is, SSD does not include any kind of<br>
> evolution in time. With this in mind, does the solver assume, when we set<br>
> the use of a mass_matrix, that it has to use optical flow instead of simply<br>
> assessing the error given by the warping? If this is true, it would be great<br>
> to see the original optimization problem, or at least the EL equations and<br>
> how they are penalized. I am very intrigued because my (very) simple python<br>
> implementation uses no time evolution and a plain fixed point scheme, and no<br>
> matter how much I beg, it does not converge to anything that lungs like a<br>
> registration. Thanks again.<br>
><br>
> Best regards<br>
><br>
> On 24 March 2017 at 15:55, Matt McCormick <<a href="mailto:matt.mccormick@kitware.com">matt.mccormick@kitware.com</a>><br>
> wrote:<br>
>><br>
>> Hi Nicolás,<br>
>><br>
>> On slide 5, E_D and E_S are referring to the optical flow registration.<br>
>><br>
>> On slide 6 is an objective function for variational registration, and<br>
>> FEM registration uses the integral of the deformation (linear elastic<br>
>> internal strain energy) as a regularization term.<br>
>><br>
>> The similarity term can vary:<br>
>><br>
>><br>
>> <a href="https://github.com/InsightSoftwareConsortium/ITK/blob/aecfac233e8815cdd0121fd2351dd3fef80d2e1b/Modules/Registration/FEM/include/itkFEMRegistrationFilter.hxx#L222-L255" rel="noreferrer" target="_blank">https://github.com/<wbr>InsightSoftwareConsortium/ITK/<wbr>blob/<wbr>aecfac233e8815cdd0121fd2351dd3<wbr>fef80d2e1b/Modules/<wbr>Registration/FEM/include/<wbr>itkFEMRegistrationFilter.hxx#<wbr>L222-L255</a><br>
>><br>
>> HTH,<br>
>> Matt<br>
>><br>
>> On Fri, Mar 24, 2017 at 11:39 AM, Nicolás Barnafi <<a href="mailto:nabw91@gmail.com">nabw91@gmail.com</a>><br>
>> wrote:<br>
>> > Hi Matt,<br>
>> ><br>
>> > Thanks the quick answer and welcome :) . That is indeed the ppt. There<br>
>> > they<br>
>> > define two energies and then write an integral in terms of some<br>
>> > "similarity"<br>
>> > which wasn't really defined. Maybe it is the E_D term? Hopefully the<br>
>> > same<br>
>> > applies for smoothness (E_S), but I'm not sure. In the code I had<br>
>> > trouble<br>
>> > because RunRegistration uses MultiResSolve, and there they instantiate a<br>
>> > solver which uses Update(). I couldn't find Update in either .h or .hxx<br>
>> > files, and apparently update is redirected (somehow) to AssembleKandM.<br>
>> > The<br>
>> > Doxygen page for the CrankNicolson solver shows a formulation of the<br>
>> > method<br>
>> > but not for what problem, and that leaves other questions without<br>
>> > answer,<br>
>> > such as the boundary conditions, elements used for the FEM formulation,<br>
>> > etc.<br>
>> > Thanks so much for your time.<br>
>> ><br>
>> > Best regards<br>
>> ><br>
>> > On 23 March 2017 at 18:34, Matt McCormick <<a href="mailto:matt.mccormick@kitware.com">matt.mccormick@kitware.com</a>><br>
>> > wrote:<br>
>> >><br>
>> >> Hi Nicolás,<br>
>> >><br>
>> >> Welcome to ITK!<br>
>> >><br>
>> >><br>
>> >> Is this the PowerPoint?<br>
>> >><br>
>> >> <a href="https://itk.org/CourseWare/Training/NonRigidRegistrationMethods" rel="noreferrer" target="_blank">https://itk.org/CourseWare/<wbr>Training/<wbr>NonRigidRegistrationMethods</a><br>
>> >><br>
>> >> Note that DeformableRegistration1.cxx (FEM registration) is discussed<br>
>> >> after the optical flow example.<br>
>> >><br>
>> >><br>
>> >> Some more information on the FEMRegistrationFilter can be found in its<br>
>> >> class documentation:<br>
>> >><br>
>> >><br>
>> >><br>
>> >> <a href="https://itk.org/Doxygen/html/classitk_1_1fem_1_1FEMRegistrationFilter.html" rel="noreferrer" target="_blank">https://itk.org/Doxygen/html/<wbr>classitk_1_1fem_1_<wbr>1FEMRegistrationFilter.html</a><br>
>> >><br>
>> >> Multiple metrics are available:<br>
>> >><br>
>> >><br>
>> >><br>
>> >> <a href="https://itk.org/Doxygen/html/classitk_1_1fem_1_1FEMRegistrationFilter.html#acc636f4752e592cd780503a5fbfeba82" rel="noreferrer" target="_blank">https://itk.org/Doxygen/html/<wbr>classitk_1_1fem_1_<wbr>1FEMRegistrationFilter.html#<wbr>acc636f4752e592cd780503a5fbfeb<wbr>a82</a><br>
>> >><br>
>> >><br>
>> >> In general, "the code reveals all" and should be given the highest<br>
>> >> amount of trust.<br>
>> >><br>
>> >><br>
>> >> Hope this helps,<br>
>> >> Matt<br>
>> >><br>
>> >> On Thu, Mar 23, 2017 at 4:32 PM, Nicolás Barnafi <<a href="mailto:nabw91@gmail.com">nabw91@gmail.com</a>><br>
>> >> wrote:<br>
>> >> > Hi everyone,<br>
>> >> ><br>
>> >> > I have been trying to understand exactly what is happening in the<br>
>> >> > DeformableRegistration1.h file without much success. The problem is:<br>
>> >> > according to the software guide, the problem being solved comes from<br>
>> >> > the<br>
>> >> > variational problem given by<br>
>> >> ><br>
>> >> > min D[image1, image2; u] + S[u]<br>
>> >> ><br>
>> >> > where D is just the SSD (or the L2 norm of (Im1 - Im2 o phi), with<br>
>> >> > phi<br>
>> >> > the<br>
>> >> > unknown displacement field) and S is a linear elastic potential. From<br>
>> >> > here<br>
>> >> > you get the euler lagrange equations (asuming some unspecified<br>
>> >> > boundary<br>
>> >> > ocndition) and solve it using some semi implicit newton-raphson<br>
>> >> > scheme.<br>
>> >> > This<br>
>> >> > is where it starts getting blurry, because the ITK ppt on deformable<br>
>> >> > registration first shows an optical flow formulation, which would<br>
>> >> > mean<br>
>> >> > that<br>
>> >> > the SSD metric isn't really what is being used, and also if I dig<br>
>> >> > deeper<br>
>> >> > in<br>
>> >> > the code, I find actually a Crank-Nicolson scheme being used, which<br>
>> >> > really<br>
>> >> > implies some kind of temporality that really does not exist in the<br>
>> >> > variational formulation. The only hint I have found was in<br>
>> >> > Modersitzki's<br>
>> >> > book where a fixed point scheme is artificially stabilized:<br>
>> >> ><br>
>> >> > A(u[k+1]) = f_u[k]<br>
>> >> ><br>
>> >> > => u[k+1] + t A(u[k+1]) = t f_u[k] + u[k].<br>
>> >> ><br>
>> >> > I would want to know what is exactly happening in that example to be<br>
>> >> > able to<br>
>> >> > validate an example I implemented in python. Thanks for your time.<br>
>> >> ><br>
>> >> > Best regards<br>
>> >> ><br>
>> >> ><br>
>> >> ><br>
>> >> > --<br>
>> >> > Nicolás Alejandro Barnafi Wittwer<br>
>> >> ><br>
>> >> > ______________________________<wbr>_________________<br>
>> >> > Community mailing list<br>
>> >> > <a href="mailto:Community@itk.org">Community@itk.org</a><br>
>> >> > <a href="http://public.kitware.com/mailman/listinfo/community" rel="noreferrer" target="_blank">http://public.kitware.com/<wbr>mailman/listinfo/community</a><br>
>> >> ><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> > --<br>
>> > Nicolás Alejandro Barnafi Wittwer<br>
><br>
><br>
><br>
><br>
> --<br>
> Nicolás Alejandro Barnafi Wittwer<br>
______________________________<wbr>_________________<br>
Community mailing list<br>
<a href="mailto:Community@itk.org">Community@itk.org</a><br>
<a href="http://public.kitware.com/mailman/listinfo/community" rel="noreferrer" target="_blank">http://public.kitware.com/<wbr>mailman/listinfo/community</a><br>
</div></div></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div><div dir="ltr">Nicolás Alejandro Barnafi Wittwer</div></div></div></div>
</div>