This project seeks to apply curvature estimation techniques to real-world scanned 3D surfaces to understand their shape and what significance the curvatures may have in 3D shape understanding.

The curvature of one point on a surface is decided by both the position of this point and curvature direction. For each point on a given surface, I need to figure out both the maximum curvature value with its direction and the minimum curvature value with its direction. I use symbol Kmax and Kmin to express the curvature and Vmax and Vmin to express the direction vector. Since we don't know the function of the surface and the points on the surface are not continuous, we must use some method for discrete math.

The first method I use is called Fitting Quadric Estimation. This method assume that the surface can be expressed as a quadric equation. If the parameters of this equation can be figured out, the K and V can be expressed by these parameters. This algorithm goes as follow:

(1) Estimate the normal vector n of the vertex and build the local principal frame with r3=n, r1=(project axis X to tangent plane), r2=r3Xr1;

(2) Map the global data to the local principal frame.

(3) Since MQ=Z, M=(x_{i}*x_{i}, x_{i}*y_{i}, y_{i}*y_{i}), Q=(a,b,c), Z=(zi), Q=(M^{t}M)^{-1} M^{t}Z. Then we get the value pf parameters a, b,c.

(4) Kmax=a+c+sqrt((a-c)^{2}+b^{2}) , Kmax=a+c-sqrt((a-c)^{2}+b^{2}), a=0.5*atan(b,a-c)

(5) local [Vmax, Vmin, Normal] = [cosa, sina, 0; -sina, cosa, 0; 0,0, 1] [r1, r2, r3]]

reverse map local [Vmax, Vmin, Normal] to global [Vmax, Vmin, Normal]

One improvement is increase to 5 parameters. [a, b, c, d, e] and z=ax^{2}+byx+cy^{2}+dx+ey while other parts remain.

The second method I use is called Laurent Series Estimation. Some steps are the same as the first one, such as estimate the normal vector, and these steps are omitted here.

First we find the curvature for each neighbor direction. Since any point on the regular surface can be expressed as Laurent series expanssion x(s)=x(0)+x'(0)s+0.5x''(0)s^{2}+O(s^{3}), and we can simplify this equation by some differential geometry method, finally we get Kp(T)=2N^{t}(q-p)/||q-p||^{2}. Here p means the point whose curvature is wanted, q means its neighbor points, N is normal vector and Kp is the curvature for vector q direction.

Then we build a 3x3 matrix M=SUM(wkTT^{t}), w is the weight of each neighbor, k is the curvature, T is the tangent projection of the vector pj-pi. The eigenvalue are [Kmax, Kmin, 0] or [Kmin, Kmax, 0], and the eigenvector matrix is [Vmax, Vmin,V] or [Vmin, Vmax, V], we don't care V.

We can see the above two method are only related with 1st order neighbors, for some reason, this could make the estimation result not reliable. So I improved the method of searching neighbors for any reasonable orders. In this experiment I used 7th order neighbor, then I combine this to method one and two as method three and four. This result may look more "smooth", however the original data are not changed.

To express the curvature for each points on surface, we paint the points and their and their neighbor mesh in different colors. But first of all, we must re-scale the curvature, since the range of the value of the curvature are very large and most values are concentrated in the lower level. We only take the values in 3 times of standard deviation for consideration and scale them from -1 to 1 by scale factor r=0.2. The RGB values are (0, 1+d, -d) for negative curvature and (d, 1-d, 0) for positive ones, here d is the scaled curvature.

The experiment result for three objects are showed as follow: