[vtkusers] Computing curvature of a surface

[Philipp.Batchelor]

Yes,

curvatures are extremely sensitive to triangulations, in particular the simple definition
that I use, but in any case, as a second order measure of the surface, it is dependent
on small bumps etc. It's one of the reasons I'd like to add some of the other more sophisticated
computations from the literature, but even then, the problem is more fundamental than which formula.
(Mathematically, a surface, and its smoothing might be close in C^0 topology, but very far in C^2
topology.)
An extreme example: an icosahedron has by definition the same curvatures everywhere
http://www-ipg.umds.ac.uk/p.batchelor/curvatures/icosa_H.jpg
but using the subdivision filter on the facets
http://www-ipg.umds.ac.uk/p.batchelor/curvatures/icosa_refined_H.jpg

We had similar problems (see the link to Andy's PhD), some suggestions based on past experience:
Basically, it requires a quite heavy postprocessing: of course, this is not really the answer to
your question
, as you want to compare the surface before and after postprocessing. It's interesting that the
decimation preserves maxima and minima.
--a very good smoothing is important for the curvatures to make sense. We used the Surface Evolver,
using it also to interactively change the triangulation until the number of triangles was roughly
what we wanted, the Evolver surface relaxation, and some other things available in Evolver. 
--the Gauss Curvature scales like [1/m^2]=1/area, the mean like [1/m]=1/distance, so a global
scaling using surface areas (scaling chosen so that they have same surface area) can be useful.
--even after smoothing, the curvatures often still have this 'leopard' look. A good improvement came
by averaging the curvatures locally, by which I mean: 
        curvature at vertex = average of curvatures at neighbouring vertex
repeated 2 or 3, up to 5 times. This widens the range over which the curvature is computed, and
makes quite good improvement. We decided on this as a simple way to tell the curvature to use 'next
to next' vertices for example.  A torus example:
http://www-ipg.umds.ac.uk/p.batchelor/curvatures/torus_notsmoothed_H.jpg
http://www-ipg.umds.ac.uk/p.batchelor/curvatures/torus_smoothed5_H.jpg
The visual appearance is more convincing, but this is sadly not the end of the story...


Ph

Jan Ehrhardt wrote:
> 
> Moin Moin,
> 
> thanks very much for any help. I have tried Phillipps classes with a
> pre-4.0 version
> and it works fine. Special thanks to Phillipp !
> 
> Now I have a new problem:
> I have a surface-representation of an object, i.e. generated with
> marching cubes and a
> second surface-representation after applying the triangle decimation
> algorithm.
> Now I compute the mean curvatures for both surfaces and the curvatures
> significantly differ
> (but the  minimal and maximal mean curvature is the same for both
> surfaces).
> 
> In http://www.medinf.mu-luebeck.de/~ehrhardt/images/vhf_curvature_dc.jpg
> and
> http://www.medinf.mu-luebeck.de/~ehrhardt/images/vhf_curvature_full.jpg
> a visualization of the curvatures with scalar-range [-0.1, 0.1] is
> schown.
> 
> Is it possible to compute something like a "normalized curvature"
> independent from the
> triangle resolution and  small local variations (maybe caused by
> discretisation artefakts).
> 
> Regards,
>         Jan
> 
> "Philipp.Batchelor" wrote:
> >
> > Yes,
> >
> > it's a funny coincidence that I hadn't been touched my curvature
> > code for a long time, but have been working this week on  the  changes to make it compatible
> > with the way vtk4.x handles Scalars. (Not without struggle, as my post
> > on Monday shows). Also, I'm trying to add a new method to compute curvatures, based
> > on principal curvatures, which also computes a curvature tensor, but this is untested yet.
> > Note that the filter was written quite a while ago, when I was still learning
> > vtk, and is certainly not up the vtk coding standards.
> >
> > In short: Gauss curvature is the angle defect at a vertex.
> >           Mean curvature, edge length x dihedral angle at edge.
> >
> > You can have a look at the main loops is the method GetK() and GetH()
> > of the filter, and copy-paste that if it is more convenient.
> >
> > ----------------------------------
> > Dr Philippe Batchelor
> > King's College London
> > London  SE1 9RT
> >
> > phone:  +44 207 955 4223
> > fax:    +44 207 955 4532
> > ----------------------------------
> >
> > "Andrew J. P. Maclean" wrote:
> > >
> > > It has all been done for you!
> > >
> > > Look at:
> > > http://www-ipg.umds.ac.uk/p.batchelor/curvatures/curvatures.html
> > >
> > > Philipp Batchelor also supplies some vtk Classes to do it:
> > > His comments are:
> > > " This a VTK C++ class, a subclass of vtkPolyDataToPolyDataFilter, as
> > > input it takes a triangulated PolyData, and the output is the same
> > > PolyData, with the desired curvatures a scalars. This is chosen by the
> > > method vtkCurvatures::SetCurvatureType(i), with i = 0 for , i = 1 for
> > > for . i = 2 should give both, but I haven't used it very much in that
> > > way. Note thata as the computation of angles uses the acos, asin, and
> > > atan2 functions, the curvatures are double.
> > > Warning As mentioned above: the mean curvature requires an oriented
> > > surface, so typically, the surface should have gone through a
> > > vtkTriangleFilter and vtkPolyDataNormals. This is quite restrictive, in
> > > particular none of the definitions really requires triangles, I just did
> > > it for convenience. If you'd want to suggest modifications, don't
> > > hesitate to email me.
> > >
> > > I use the class in a local library, so that the curvature class can be
> > > used in Tcl. I normally use VTK 3.1, on a Sun Ultra-30, with the Sun
> > > compiler.
> > >
> > > As they are relatively small files, they are provided directly: Download
> > > vtkCurvatures.h and DownloadvtkCurvatures.cxx "
> > >
> > > Andrew
> > >
> > > ___________________________________________
> > > Andrew J. P. Maclean
> > > Postal:
> > > Australian Centre for Field Robotics
> > > The Rose Street Building J04
> > > The University of Sydney  2006  NSW
> > > AUSTRALIA
> > >
> > > Room:  106
> > > Phone: +61 2 9351 3283
> > > Fax:   +61 2 9351 7474
> > >        http://www.acfr.usyd.edu.au/
> > >
> > > ___________________________________________
> > >
> > > -----Original Message-----
> > > From: vtkusers-admin at public.kitware.com
> > > [mailto:vtkusers-admin at public.kitware.com] On Behalf Of Jan Ehrhardt
> > > Sent: Wednesday, 3 July 2002 18:11
> > > To: vtkusers
> > > Subject: [vtkusers] Computing curvature of a surface
> > >
> > > hi folks,
> > >
> > > I have the following problem: I want to register two polygonal surfaces
> > > using nonlinear transformations.
> > > I have implemented an algorithm similar to andresens and nielsen
> > > "geometry constrained diffusion".
> > > This works like ICP but for nonlinear transformations.
> > > Actually, I match the nearest points on the two surfaces.
> > > The idea is not only to use the spatial distance but also the normal
> > > vektors and curvature of the points
> > > to identify corresponding points.
> > > The normals caan I compute using vtkPolyDataNormals, but how to compute
> > > the curvature of the surface points?
> > > Concrete:
> > > I have an polydata surface (no associated image or scalar values!).
> > > How can I compute the curvature for each point (triangle vertex) from
> > > the polydata ?
> > > Curvature means min and max curvature and /or gaussian curvature.
> > >
> > > Regards and many thanks for ANY help,
> > >
> > > Jan
> > >
> > > _______________________________________________
> > > This is the private VTK discussion list.
> > > Please keep messages on-topic. Check the FAQ at: <http://public.kitware.com/cgi-bin/vtkfaq>
> > > Follow this link to subscribe/unsubscribe:
> > > http://public.kitware.com/mailman/listinfo/vtkusers
> >
> > --

-- 
----------------------------------
Dr Philippe Batchelor
King's College London  	
London	SE1 9RT

phone: 	+44 207 955 4223
fax:	+44 207 955 4532
------------------------------------------------------