# geometry-processing-curvature
**Repository Path**: mirrors_lepy/geometry-processing-curvature
## Basic Information
- **Project Name**: geometry-processing-curvature
- **Description**: No description available
- **Primary Language**: Unknown
- **License**: MPL-2.0
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No
## Statistics
- **Stars**: 0
- **Forks**: 0
- **Created**: 2021-10-22
- **Last Updated**: 2023-08-19
## Categories & Tags
**Categories**: Uncategorized
**Tags**: None
## README
# Geometry Processing - Curvature
> **To get started:** Clone this repository then issue
>
> git clone --recursive http://github.com/[username]/geometry-processing-curvature.git
>
## Installation, Layout, and Compilation
See
[introduction](http://github.com/alecjacobson/geometry-processing-introduction).
## Execution
Once built, you can execute the assignment from inside the `build/` by running
on a given mesh:
./curvature [path to mesh.obj]
## Background
In this assignment we explore discrete curvature quantities computed on a
surface. These quantities give us local information about a shape. Beyond
inspecting the surface (the extent of this assignment), these quantities become
the building blocks to:
- define energies to minimize during smoothing/deformation,
- identify salient points and curves on the shape, and
- provide initial conditions/constraints for _remeshing_.
The fundamental difference between a segment on the real line and a curve is
the introduction of [curvature](https://en.wikipedia.org/wiki/Curvature). This
is quite natural and intuitive. When we draw a 1D object in the plane or in
space we have the freedom to let that object bend. We quantify this "bending"
locally as curvature.
[Curvature](https://en.wikipedia.org/wiki/Curvature#Surfaces) is also the
fundamental difference between a chunk (i.e., subregion) of the [Euclidean
Plane](https://en.wikipedia.org/wiki/Plane_(geometry)) and a
[surface](https://en.wikipedia.org/wiki/Surface_(mathematics)) that has been
[immersed](https://en.wikipedia.org/wiki/Immersion_(mathematics)) in 
(or
elsewhere). Unlike curves, surfaces can bend in each direction at any point.
We start our discussion assuming a smooth surface 
. We would like to
categorize points on the surface 
in terms of how the surface bends or
curves locally.
### Curvature of planar curves
Let us briefly recall how
[curvature](https://en.wikipedia.org/wiki/Curvature#Precise_definition) is
defined for a [planar curve](https://en.wikipedia.org/wiki/Plane_curve)
.
There are multiple equivalent definitions.
#### Osculating circle
We can define the [tangent](https://en.wikipedia.org/wiki/Tangent) direction at
a point 
as the limit of the
[secant](https://en.wikipedia.org/wiki/Secant_line) formed between 
and
another point on the curve 
as 
approaches :
It always possible, and often convenient, to assume without loss of generality
that 
is an [arc length
parameterization](https://en.wikipedia.org/wiki/Arc_length) of the curve 
so
that 
and therefor the unit tangent vector is simply 
.
In an analogous fashion, we can consider the limit of the
[circumcircle](https://en.wikipedia.org/wiki/Circumscribed_circle)
that passes
through and points 
and 
before and after it on the curve:
This limit circle is called the [osculating
circle](https://en.wikipedia.org/wiki/Osculating_circle) at the point on
the curve . By construction the tangent of the curve and the circle match at
: they're both 
. The
[radius](https://en.wikipedia.org/wiki/Radius_of_curvature) 
of the
osculating circle 
at the the point is proportional to how straight
the curve is locally: as the curve becomes more and more straight then the
radius tends toward infinity. This implies that the radius is inversely
proportional to the "curvy-ness" of the curve. Hence, the inverse of the radius
is dubbed the curvature:
The radius is a non-negative measure of length with units meters, so the
curvature 
is an non-negative scalar with units 1/meters. The radius of the
osculating circle can also be written as a limit of the [circumcircle
radius](https://en.wikipedia.org/wiki/Circumscribed_circle#Cartesian_coordinates_from_cross-_and_dot-products):
#### Signed curvature
Plugging in our arc-length parameterization this reveals that the curvature
(inverse of radius) is equal to the magnitude of change in the tangent or
equivalently the magnitude of second derivative of the curve:
Because we chose the arc-length parameterization, the only change to the
tangent vector is a change in _direction_ (as opposed to magnitude, since
). This means that the change--as a vector itself--is
_orthogonal_ to the tangent. In other words, the change in tangent 
points
along the normal direction 
:
If we define an orientation to our curve then we can endow the curvature with
a [sign](https://en.wikipedia.org/wiki/Sign_(mathematics)) based on whether the
center of the osculating circle lies on the [left or right
side](https://en.wikipedia.org/wiki/Right-hand_rule) of the curve. As already
established, the tangent of the osculating circle and the curve agree, so the
vector pointing toward the circle's center must be
[perpendicular](https://en.wikipedia.org/wiki/Perpendicular) to the tangent:
i.e., in either the positive or negative
[normal](https://en.wikipedia.org/wiki/Normal_(geometry)) directions.
If the orientation agrees with increasing the arc-length parameter , then the sign can
be
determined by comparing the second derivative vector to the unit normal
. The [_**signed
curvature**_](https://en.wikipedia.org/wiki/Curvature#Signed_curvature) at a point is thus given by:
#### Moving point analogy
This definition neatly conforms to our intuition of a curve as the trajectory
of a moving point. Imagine the curved formed by driving along a particular
trajectory 
, where we really interpret 
as time.
While 
corresponds to your velocity vector and 
corresponds to
your speed, the arc-length (re-)parameterization would correspond to having
your friend re-trace your path traveling at a perfectly uniform speed 
, where your friends "time" may be different from yours (it may take
longer or shorter depending if you drove fast or slow).
Curvature in the path corresponds to _turning_ and quite literally the amount
by which your friend needs to turn the steering wheel away from the "straight"
position: on a straight course, the steering wheel remains at zero-angle
position and the curvature is zero, on a circular course the steering wheel is
fixed at a constant angle in the left or right direction corresponding to
constant positive or negative curvature respectively.
Changing the steering wheel changes the _direction_ of the vehicle's velocity.
For your friend driving at constant speed, this is the _only_ change admissible
to the velocity, hence the curvature exactly corresponds to 
and to the
steering wheel angle.
> If somebody wants to make a Sega [Out
> Run](https://en.wikipedia.org/wiki/Out_Run) inspired gif showing a steering
> wheel turning next to a little car tracing a curve, I'll be very impressed.
#### Turning number
The integrated signed curvature around a [closed
curve](https://en.wikipedia.org/wiki/Curve) must be an integer multiple of
:
where 
is an integer called the "turning number" of the curve.
This is a bit surprising at first glance. However, in the _moving point
analogy_ a closed curve corresponds to a period trajectory (e.g., driving
around a race-track). When we've made it once around the track, our velocity
direction (e.g., the direction the vehicle is facing) must be pointing in the
original direction. That is, during the course, the car either have turned all the
way around once (
) or turned as much clockwise and it did
counter-clockwise (e.g., on a figure 8 course: 
), or made multiple
loops, etc.
#### Discrete curvature
In the discrete world, if a curve is represented as a piecewise-linear chain of
segments, then it's natural to associate curvature with vertices: the segments
are flat and therefor contain no curvature.
A natural analog to the definition of curvature as
the derivative of the tangent vector
(i.e., 
) is to define _discrete curvature_ as the change in
tangent direction between discrete segments meeting at a vertex:
that is, the signed [_exterior
angle_](https://en.wikipedia.org/wiki/Internal_and_external_angles) 
at
the vertex 
.
The turning number theorem for continuous curves finds an _immediate_ analog in
the discrete case. For a closed polygon the discrete signed angles must sum up
to a multiple of in order to close up:
In this way, we _preserve the structure_ found in the continuous case in our
discrete analog. This structure preservation leads to an understanding of the
exterior angle as an approximation or discrete analog of the _locally
integrated_ curvature.
Alternatively, we could literally fit an circle to the discrete curve based on
local samples and approximate curvature as the inverse radius of the osculating
circle. This curvature measure (in general) will not obey the turning number
theorem, but (conducted properly) it will converge to the pointwise continuous
values under refinement (e.g., as segment length shrinks).
We will explore these two concepts for surfaces, too: discrete analogs that
preserve continuous structures and discretizations that approximate continuous
quantities in the limit.
### Curvature(s) on surfaces
A surface can be curved locally in multiple ways. Consider the difference
between a flat piece of paper, a spherical ping-pong ball and a saddle-shaped
[Pringles chip](https://en.wikipedia.org/wiki/Pringles). The Pringles chip is
the most interesting because it curves "outward" in one direction and "inward"
in another direction. In this section, we will learn to distinguish and
classify points on a surface based on how it curves in each direction.
#### Normal curvature
The simplest way to extend the curvature that we defined for planar curves to a
surface is to _slice_ the surface through a given point 
with a
[plane](https://en.wikipedia.org/wiki/Plane_(geometry)) 
that is parallel
to the [surface normal](https://en.wikipedia.org/wiki/Normal_(geometry))
.
The (local) intersection of the surface and the plane will trace a
curve , upon which we can immediately use the planar curvature definition
above.
](images/normal-curvature.svg)
There are infinitely many planes that pass through a given point and lie
parallel to a given normal vector : the plane can rotate around the
normal by any angle 
. For each choice of , the plane will define
an intersecting curve 
and thus for every angle there will be a
[_normal curvature_](https://en.wikipedia.org/wiki/Curvature#Normal_sections):
#### Mean curvature
Normal curvature requires choosing an angle, so it doesn't satiate our
desire to reduce the "curvy-ness" to a single number for any point on the
surface. A simple way to reduce this space of normal curvatures is to, well,
average all possible normal curvatures. This defines the [mean
curvature](https://en.wikipedia.org/wiki/Mean_curvature):
#### Maximum and minimum curvature
Another obvious way to reduce the space of normal curvatures to a single number
is to consider the maximum or minimum normal curvature over all choices of :
Collectively, these are referred to as the [principal
curvatures](https://en.wikipedia.org/wiki/Principal_curvature) and
correspondingly the angles that maximize and minimize curvature are referred to
as the principal curvature directions:
[Euler's
theorem](https://en.wikipedia.org/wiki/Euler%27s_theorem_(differential_geometry))
states that the normal curvature is a quite simple function of and the
principal curvatures:
([proof](https://math.stackexchange.com/a/1783316/35376)).
There are two immediate and important consequences:
1. the principal curvature directions (
and 
) are orthogonal, and
2. the mean curvature reduces to the average of principal curvatures:
> For more theory and a proof of Euler's theorem, I recommend "Elementary
> Differential Geometry" by Barret O'Neill, Chapter 5.2.
#### Gaussian curvature
Maximum, minimum and mean curvature reduce curvature to a single number, but
still cannot (alone) distinguish between points lying on a round ping-pong
ball, a flat sheet of paper, the cylindrical Pringles can and a saddle-shaped
Pringles chip.
The neck of this cartoon elephant--like a Pringles chip--bends inward in one
direction (positive 
) and outward in the other
direction (negative 
).
Figure Caption: Maximum 
, minimum 
, and Gaussian curvature 
.
The _product_ of the principal curvatures maintains the disagreement in sign
that categories this saddle-like behavior. This product is called [Gaussian
curvature](https://en.wikipedia.org/wiki/Gaussian_curvature):
#### Relationship to surface area
Both mean and Gaussian curvature have meaningful relationships to surface
area.
##### Mean Curvature as area gradient
Let us consider a seemingly unrelated yet familiar problem. Suppose we would
like to _flow_ a given surface in the direction that shrinks its surface area.
That is, we would like to move each surface point in the direction that
minimizes surface area.
The surface area of may be written as an integral of unit density:
There are many expressions that 
. We can choose an expression that is
especially easy to work with. Namely, the small change in position over a small
change in position is a unit vector.
The norm of the gradient is a non-linear function involving square roots, but
since the magnitude is one then the squared magnitude is also one (
. This allows us to write the surface area as a quadratic function of
positions and familiarly as the Dirichlet energy:
By abuse of notation we can say that 
is a functional (function
that takes a function as input) and measures the surface area of the surface
defined by the embedding function 
. Now, let's consider the
[functional derivative](https://en.wikipedia.org/wiki/Functional_derivative) of
with respect to . This special type of derivative can be written
as:
where 
is an _arbitrary_ function. That is, we consider the limit of
a tiny perturbation of the function in any way.
We can identify this limit by considering the derivative of the perturbation
magnitude 
evaluated at zero:
Feeding in our Dirichlet energy definition of we can start
working through this derivative:
Assuming that is closed (no boundary), then applying [Green's identity](https://en.wikipedia.org/wiki/Green%27s_identities#Green's_first_identity) leaves us with:
This still leaves us with an expression of the derivative written as an integral
involving this arbitrary function . We would like to have a more
compact expression to evaluate 
at some query point
on the surface.
Since this must be true for any choice of perturbation function , we
can choose to be a function that is 
everywhere on the domain
except in the region just around 
, where makes a little
"bump" maxing out at 
. Since this bump can be made
arbitrarily skinny, we can argue that can be factored out of the
integral above (if 
everywhere except 
arbitrarily close to
, then the integral just evaluates to 
at
):
This reveals to us that the Laplacian of the embedding function indicates the
direction and amount that the surface should move to decrease surface area.
The Laplacian 
of a function 
on the surface does not depend on the
choice of parameterization. It is defined as the divergence of the gradient of
the function or equivalently the trace of the Hessian:
If we generously choose 
and 
to vary in the principal directions
and above. In this case, the Laplacian of the position function
reduces to the sum of principal curvatures times the normal (recall the
definition of [curvature normal](#curvature-normal)):
where 
is called the _**mean curvature normal**_ vector. We have
shown that the mean curvature normal is equal half the Laplacian of the
embedding function, which is in turn the gradient of surface area.
##### Gaussian Curvature as area distortion
As the product of principal curvatures, Gaussian curvature measures zero
anytime one (or both) of the principal curvatures are zero. Intuitively, this
happens only for surfaces that curve or bend in one direction. Imagine rolling
up a sheet of paper. Surfaces with zero Gaussian curvature 
are called
_developable surfaces_ because the can be flattened (developed) on to the flat
plane (just as you might unroll the piece of paper) _without_ stretching or
shearing. As a corollary, surfaces with non-zero Gaussian curvature _cannot_ be
flattened to the plane without stretching some part.
Locally, Gaussian curvature measures how far from developable the surface is:
how much would the local area need to stretch to become flat.
First, we introduce the [Gauss map](https://en.wikipedia.org/wiki/Gauss_map), a
continuous map 
from every point on the surface to the unit
sphere 
so that 
, the unit normal at .
Consider a small patch on a curved surface. Gaussian curvature 
can
equivalently be defined as the limit of the ratio between the area
area _swept_ out by the unit normal on the Gauss map 
and
the area of the surface patch :
Let's consider different types of regions:
- flat: 
because the Gauss map is a point,
- cylindrical: because the Gauss map is a curve,
- spherical: 
because the Gauss map will maintain positive swept-area,
and
- saddle-shaped: 
because the area on the Gauss map will maintain
_oppositely_ oriented area (i.e., from the spherical case).
Similar to the turning number theorem for curves, there exists an analogous
[theorem for surfaces](https://en.wikipedia.org/wiki/Gauss-Bonnet_theorem)
stating that the total Gaussian curvature must be an integer multiple of :
where 
is the [Euler
characteristic](https://en.wikipedia.org/wiki/Euler_characteristic) of the
surfaces (a topological _invariant_ of the surface revealing how many
[holes](https://en.wikipedia.org/wiki/Genus_(mathematics)) the surface has).
In stark contrast to mean curvature, this theorem tells us that we cannot add
Gaussian curvature to a surface without:
1. removing an equal amount some place else, or
2. changing the topology of the surface.
Since changing the topology of the surface would require a discontinuous
deformation, adding and removing Gaussian curvature must also balance out for
smooth deformations. This simultaneously explains why a cloth must have
wrinkles when draping over a table, and why a deflated basketball will not lie
flat on the ground.
#### Shape operator
There is yet another way to arrive at principal, mean and Gaussian curvatures.
Consider a point on a surface with unit normal vector 
. If we
pick a unit tangent vector 
(i.e., so that 
), then we can ask
how does the normal change as we move in the direction of along the
surface:
we call 
the [_**shape
operator**_](https://en.wikipedia.org/wiki/Differential_geometry_of_surfaces#Shape_operator)
at the point . Just as how in the definition of [curvature normal](#curvature-normal), the
curvature normal must point in the normal direction, the shape operator takes
as input a tangent vector and outputs another tangent vector (i.e., the change
in the unit normal must be tangential to the surface; no change can occur in
the normal direction itself).
Locally, the tangent vector space is two-dimensional spanned by basis vectors
so we can think of the
shape operator as a mapping from 
to . As a differential operator,
the shape operator is a _linear operator_. This means we can represent its
action on a tangent vector 
as a matrix:
Given 
and 
are the principal curvature directions (as unit 2D tangent
vectors) we can rotate our coordinate frame to align 
and 
with the
principal curvature directions. The shape operator takes on a very special
form:
> Consider why the off-diagonal terms are zero. Think about the _extremality_
> of the principal curvatures.
We have actually conducted an [eigen
decomposition](https://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix) on
the shape operator. Reading this progression backwards, the eigen decomposition
of the shape operator expressed in any basis will reveal:
1. the principal curvatures as the eigen values, and
2. the principal curvature directions as the eigen vectors.
### Discrete curvatures on surfaces
#### Discrete mean curvature normal via discrete Laplace
By now we are very familiar with the discrete Laplacian for triangle meshes:
where 
are the mass and cotangent matrices respectively.
When applied to the vertex positions, this operator gives a point-wise (or
rather integral average) approximation of the mean curvature normal:
Stripping the magnitude off the rows of the resulting matrix would give the
_unsigned_ mean curvature. To make sure that the sign is preserved we can check
whether each row in 
agrees or disagrees with consistently oriented
per-vertex normals in 
.
This connection between the Laplace operator and the mean curvature normal
provides additional understanding for its use as a geometric smoothing operator
(see "Computing Discrete Minimal Surfaces and Their Conjugates" [Pinkall and
Polthier 1993]).
#### Discrete Gaussian curvature via angle defect
On a discrete surface represented as a triangle mesh, curvature certainly can't
live on the flat faces. Moreover, Gaussian curvature can't live along edges
because we can always _develop_ the triangles on either side of an edge to the
plane without stretching them. In fact we can develop any arbitrarily long
chain of faces connected by edges so long as it doesn't form a loop or contain
all faces incident on a vertex. This hints that discrete Gaussian curvature
(like curvature for curves) must live at vertices.
Using the definition of Gaussian curvature in terms of the area on the [Gauss
map](#gauss-map), flat faces correspond
points on the Gauss map (contributing nothing), edges correspond to area-less
curves (traced by their [dihedral
angles](https://en.wikipedia.org/wiki/Dihedral_angle)), but vertices correspond
to spherical polygons connecting face normal-points. The area 
subtended on
the Gauss map is call the [solid
angle](https://en.wikipedia.org/wiki/Solid_angle). Conveniently, this area is
simply the [angle
defect](https://en.wikipedia.org/wiki/Angular_defect#Descartes.27_theorem) of
internal angles 
incident on the 
-th vertex contributed by each -th
incident face:
!["Gaussian Curvature and Shell Structures" [Calladine
1986]](images/angle-defect.png)
Thus, our discrete analog of locally _integrated_ Gaussian curvature is given
as the angle defect at the -th vertex. The local integral average (or
_pointwise_) discrete Gaussian curvature is the angle defect divided by the
local area associated with the -th vertex:
By way of closing up the Gauss map, closed polyhedral surfaces (i.e., meshes)
will obey the
[Gauss-Bonnet](https://en.wikipedia.org/wiki/Gauss-Bonnet_theorem)
[above](#gauss-bonnet), too:
We can connect this to [Euler's
formula](https://en.wikipedia.org/wiki/Euler_characteristic) for polyhedra in our very first
assignment:
where 
are the number of vertices, edges and faces respectively.
#### Approximation and eigen decomposition of the shape operator
Alternatively, we can approximate all curvatures of a surface by locally
fitting an analytic surface and _reading_ off its curvature values. Since
planes have no curvature, the simplest type of analytic surface that will give
a non-trivial curvature value is a quadratic surface.
Thus, the algorithm proceeds as follows. For each vertex of the given mesh,
1. gather a sampling of points in the vicinity. For simplicity, let's just
grab all other vertices that share an edge with or share an edge with a
vertex that shares an edge with (i.e., the "two-ring" of ). For most
sane meshes, this will provide enough points. Gather the positions of these
points _relative_ to (i.e., 
) into a matrix 
.
2. Next, we are going to define a quadratic surface as a height field above
some two-dimensional plane passing through . Ideally, the plane is
orthogonal to the normal at . To find such a plane, compute the
[principal-component
analysis](https://en.wikipedia.org/wiki/Principal_component_analysis) of
(i.e., conduct eigen decomposition on 
). Let 
be the coefficients for two most principal directions (call them the -
and - directions) corresponding to each point in , and let 
be the "height" of each point in the least principal direction (call
it the 
-direction).
3. An quadratic function as a height-field surface passing through the origin
is given by:
We have sets of 
values and values. Treat this as a
least-squares fitting problem and solve for the 5 unknown coefficients.
(`igl::pinv` is good for solving this robustly).
4. Each element of the shape operator for the graph of a quadratic function
over the plane has a closed form expression. You need to derive these by hand.
Just kidding. The shape operator can be constructed as the product of two
matrices:
known as the second and first fundamental forms respectively. The entries of
these matrices categorize the stretch and bending in each direction:
See Table 1 of "Estimating Differential Quantities Using Polynomial Fitting of
Osculating Jets" [Cazals & Pouget 2003] to double check for typos :-).
5. Eigen decomposition of 
reveals the principal curvatures and
_and_ the principal tangent directions (in the 
PCA basis).
6. Lift the principal tangent directions back to world 
coordinates.
## Tasks
[Download](https://archive.org/details/ElementaryDifferentialGeometry) Barret
O'Neill's book. This is my go-to differential geometry book. The section on
curvature and the shape operator should help resolve questions and fill in
missing proofs above.
### Blacklist
- `igl::gaussian_curvature`
- `igl::internal_angles` (or any of the other overloads)
- `igl::principal_curvatures`
### Whitelist
- `igl::adjacency_matrix.h`
- `igl::cotmatrix`
- `igl::invert_diag`
- `igl::massmatrix`
- `igl::per_vertex_normals`
- `igl::pinv`
- `igl::slice`
- `igl::sort`
- `igl::squared_edge_lengths`
### `src/mean_curvature.cpp`
Compute the discrete mean curvature at each vertex of a mesh (`V`,`F`) by
taking the signed magnitude of the mean curvature normal as a _pointwise_ (or
_integral average_) quantity.
### `src/internal_angles.cpp`
Given (squared) edge-lengths of a triangle mesh `l_sqr` compute the internal
angles at each corner (a.k.a. wedge) of the mesh.
### `src/angle_defect.cpp`
Compute the discrete angle defect at each vertex of a triangle mesh
(`V`,`F`), that is, the _locally integrated_ discrete Gaussian
curvature.
### `src/principal_curvatures.cpp`
Approximate principal curvature values and directions locally by considering
the two-ring neighborhood of each vertex in the mesh (`V`,`F`).