[Insight-developers] Laplace Beltrami filter output

Mon Feb 1 18:52:23 EST 2010

Hi Michel,

     It looks like 2010 is going to be
     the year of Numerical Libraries   :-)

There is Eigen analysis code in

      itkSymmetricEigenAnalysis.h

I'm not sure it that could manage the
type of matrices that you have.

Regarding "netlib", we have to be VERY
careful with code taken from that web site.

This code often lack proper Copyright
and License notices.

We just come from removing code that
turned out to be copyrighted by ACM
and subject to an absurdly restrictive
license.

Instead of taking code from the soup
of "netlib.org", we should go to the
original web sites.

It seems that ARPACK is at:
http://www.caam.rice.edu/software/ARPACK/

and its license seems to be a BSD one:
http://www.caam.rice.edu/software/ARPACK/RiceBSD.txt#LICENSE

So, that's an acceptable proposition.

Maybe we should add this to the discussion
about mathematical libraries in the Wiki:

http://www.itk.org/Wiki/Proposals:Sparse_Linear_Solvers

     Regards,

            Luis

------------------------------
On Mon, Feb 1, 2010 at 2:05 PM, Michael Bowers <mbowers at cis.jhu.edu> wrote:
>
>
> Arnaud,
>
>
>
> Thanks very much for checking, and thanks for the reference to the paper.  I
> am looking forward to reading it.
>
>
>
> I am sorry I’ve been incommunicado for so long.  I’ve been working on an
> issue I encountered when I implemented the filter for open surfaces.  I’m
> not sure what exactly changes about the eigensystem calculation, but
> apparently the problem is more complicated and
> vnl_sparse_symmetric_eigensystem will not produce a solution.  The error it
> generates is:
>
>
>
> Error: vnl_sparse_symmetric_eigensystem:
>
>   disastrous loss of orthogonality - internal error
>
> itk::ERROR: LaplaceBeltramiFilter(0xe85a4c0): Eigensystem failure with
> result:  -8
>
>
>
> This class uses FORTRAN routine dnslao from netlib.
>
>
>
> We have done most of our research in MATLAB, and our code there uses the
> EIGS function which calculates eigensystems using ARPACK.  I built the c++
> wrapper arpack++ and integrated it into the LaplaceBeltrami filter code, and
> I am able to produce a solution for the eigensystem that matches what we’ve
> gotten in MATLAB.  I think for this filter to work the sparse matrix
> eigenvalue solver from ARPACK may be required.
>
>
>
> Does this present a problem for development within itk?  As far as I can
> tell, the desired ARPACK routine (dsaupd) is not included in the v3p netlib
> library in vxl.  Are there other possible solutions that I’m missing?  Would
> it be possible to add the ARPACK FORTRAN code for this routine to the v3p
> netlib library?
>
>
>
> Any ideas you might have would be appreciated.
>
>
>
> Thank you,
>
>
>
> Mike B.
>
>
>
>
>
>
>
>
>
> From: Arnaud Gelas [mailto:Arnaud_Gelas at hms.harvard.edu]
> Sent: Monday, February 01, 2010 9:03 AM
> To: Michael Bowers
> Subject: Re: [Insight-developers] Laplace Beltrami filter output
>
>
>
> Hi Mike,
>
>
>
> How are you doing?
>
> I have just wanted to know how you are doing with the code?
>
>
>
> Thanks to one itk user, I have fixed one major bug on the computation of the
> mean curvature (with QuadEdgeMesh). Note that this filter implements the
> following paper http://www.geometry.caltech.edu/pubs/DMSB_III.pdf
>
> It appears to me (quite recently), that it is exactly the same formulation.
> Indeed Laplace Beltrami are related to the mean curvature.
>
> If you need some help, we can set up a web chat, web conference, to discuss
> about all of this.
>
>
>
> Arnaud
>
>
>
>
>
> On Dec 21, 2009, at 2:24 PM, Michael Bowers wrote:
>
>
>
> Arnaud,
>
>
>
> Thank you very much for the reply.  I would be grateful for any contribution
> you care to make – either comments via e-mail or directly to the code.  My
> knowledge of itk is very superficial so far so I tend to revert back to vnl
> operations when I’m sure there are more elegant itk methods.  All advice is
> greatly appreciated.
>
>
>
> I have attached the source for the filter – but are you able to check it out
> directly from the NAMIC sandbox?  It should be under
> JHU-CIS-ComputationalAnatomy.
>
>
>
> Thanks again.  Looking forward to working with you,
>
>
>
> Mike B.
>
>
>
>
>
>
>
>
>
>
>
> From: Arnaud Gelas [mailto:Arnaud_Gelas at hms.harvard.edu]
> Sent: Friday, December 18, 2009 4:35 PM
> To: Michael Bowers
> Subject: Re: [Insight-developers] Laplace Beltrami filter output
>
>
>
> Hi Michael,
>
>
>
> First, thanks for sending me the paper, I appreciate!
>
>
>
> Let me write my comments:
>
>
>
> Page 3: your expression of a(v)
>
> when the degree of a vertex becomes too large, or if you have obtuse angle,
> it is better to use more "advanced formulation". Have a look to the
> following insight journal
> paper: http://insight-journal.org/browse/publication/302. ooops, my bad you
> write later about it :-): Great!!!
>
> Note that you can find the code for computing this area in itk itself. I ll
> send you the link to the right class later.
>
>
>
> Regarding the implementation, if you are interested in I can give you some
> guidance on how to improve the code. Just let me know if it is fine for you,
> I can paste you some piece of code, and explain them to you.
>
> We can either exchange by email, or if you are interested in I can directly
> write comments or past some code in your souce. Just let me know if you feel
> comfortable with any of these ways.
>
>
>
> Arnaud
>
>
>
>
>
>
>
> On Dec 18, 2009, at 10:08 PM, Michael Bowers wrote:
>
>
>
>
> Arnaud,
>
>
>
> Thanks for your comments on my e-mail below.  I appreciate your interest in
> the code so have attached a draft of the Insight Journal Submission which
> has a detailed section describing the algorithm.  I have been reading the
> pdf file you referenced (and cited it in the paper), it is similar to what
> we’re trying to do I think.
>
>
>
> If you have any further comments or questions they would be appreciated.
>
>
>
> Thank you,
>
>
>
> Michael Bowers
>
> Center for Imaging Science
>
> Johns Hopkins University
>
>
>
>
>
>
>
> From: Arnaud GELAS [mailto:arnaud_gelas at hms.harvard.edu]
> Sent: Wednesday, December 09, 2009 6:36 PM
> To: Michael Bowers
> Cc: insight-developers at itk.org
> Subject: Re: [Insight-developers] Laplace Beltrami filter output
>
>
>
> Hi Michael,
>
> Good job!
>
> I am a bit surprised by the computation of the matrix, which is not the
> standard way. What is your reference?
>
> I am pretty sure it is connected to the expression involving the cotangent
> (see http://alice.loria.fr/publications/papers/2009/spectral_course/spectral_course.pdf for
> full reference on the topic). You may consider follow the
> what is in this reference.
>
> To fully benefit from the itk::QuadEdgeMesh, you may consider having a look
> at itk::QuadEdgeParam (some part are extremely similar, to what you want to
> do).
>
> For example you'll find some expression for each edge, $\Delta_{ij}=...$.
> You can find some functors with matrix coefficients you are interested in in
> itkQuadEdgeMeshParamMatrixCoefficients.h. If there is not exactly the same
> one, you can create your own by inheriting from itk::MatrixCoefficients (btw
> I guess it would be good to rename these classes).
>
> Instead of iterating on the face container, I would recommend to iterate on
> the edge container. Have a look at itk::QuadEdgeMesh::GetEdgeCells(), then
> you can easily iterate on all half edges.
> To get the origin of the half edge, have a look at
> itk::GeometricalQuadEdge::GetOrigin()  and GetDestination(). Since your
> matrix is symmetric, it is straightforward to fill it.
>
> Regarding now the decomposition in itself, vnl does the job, but I am not so
> sure about the performance... Check in the reference above what would be the
> best way to proceed. Especially svd methods generally compute the other way
> around than the one you are interested in (you are generally interested in
> lower frequency and you will first get the higher one). In the reference,
> they present a very efficient and nice way to do it.
>
> If you need any help, please let me know
>
> Arnaud
>
>
>
> On 12/09/2009 02:34 PM, Michael Bowers wrote:
>
>
>
> Itk Developers,
>
>
>
> I have placed in the NAMIC Sandbox (under JHU-CIS-ComputationalAnatomy) an
> initial draft of code to calculate a Laplace-Beltrami operator on the
> vertices of a QuadEdgeMesh.  The code is implemented as a
> QuadEdgeMeshToQuadEdgeMesh filter.  A QuadEdgeMesh is the input, which the
> filter copies to the output QuadEdgeMesh.  If basis functions are requested
> on the surface (by SetEigenValueCount > 0) , the filter sets the output
> Mesh’s PointData to the values of the 0th harmonic.  The user can get the
> values of the other harmonics by SetSurfaceHarmonic(N), then calling
> GetOutput, and the PointValue of the output Mesh will be set to the Nth
> surface harmonic.
>
>
>
> This behavior is illustrated in the attached test program.
>
>
>
> Luis has remarked that this is not the typical dynamic by which filters
> provide their output data.  He recommended putting the question to the itk
> developer community in general to ask for suggestions on how to modify the
> code to conform with standard itk pipeline design.
>
>
>
> Any recommendations/comments would be greatly appreciated.
>
>
>
> Thank you,
>
>
>
> Michael Bowers
>
> Center For Imaging Science
>
> Johns Hopkins University
>
>
>
>
>
> <LaplaceBeltramiFilter.pdf>
>
>
>
>