Surface Curvature from Range Images (or Depth Maps)

[Brickle Macho]

I have implemented mean and gaussian curvature methods in python and thought they may be appropriate for inclusion to skimage.  An outline of the code is provided below.  It assumes a surface defined as a function of two coordinates, e.g. z = Z(x, y).  The curvature calculations are from the following two papers

* Kurita, T., & Boulanger, P. (1992). Computation of Surface Curvature from Range Images Using Geometrically Intrinsic Weights. In MVA (pp. 389-392).

* Zhao, C., Zhao, D., & Chen, Y. (1996, August). Simplified Gaussian and mean curvatures to range image segmentation. In Pattern Recognition, 1996., Proceedings of the 13th International Conference on (Vol. 2, pp. 427-431). IEEE.

If appropriate let me know and I will read the contribution/development documentation, make the code more robust and submit formally.

Regards,

Michael.
--


Here is a basic outline of the code (H and K are used as they appear common terms used in the literature) :
---------------BEGIN CODE--------------
import numpy as np

Zy, Zx = np.gradient(Z)                                                    
Zxy, Zxx = np.gradient(Zx)                                                 
Zyy, _ = np.gradient(Zy) 

# Mean Curvature - equation (3) from Kurita and Boulanger (1992) paper
# See also Surface in 3D space, http://en.wikipedia.org/wiki/Mean_curvature
H = (1 + (Zx ** 2)) * Zyy + (1 + (Zy ** 2)) * Zxx - 2 * Zx * Zy * Zxy      
H = H / ((2 * (1 + (Zx ** 2) + (Zy ** 2))) ** 1.5)                         
                                                                               
# Gaussian Curvature - equation (4) from Kurita and Boulanger (1992) paper 
K = (Zxx * Zyy - (Zxy ** 2)) /  ((1 + (Zx ** 2) + (Zy **2)) ** 2)

# Simplified Mean Curvature - equation (3) from Zhao et.al (1996) paper               
H = Zxx + Zyy                                                              
                                                                               
#  Simplified Gaussian Curvature - equation (3) from Zhao et.al (1996) paper           
K = Zxx * Zyy - (Zxy ** 2)  
---------------END CODE--------------

[Stéfan van der Walt]

Dear Michael

On Wed, Nov 13, 2013 at 3:19 PM, Brickle Macho <brickl...@gmail.com> wrote:
> I have implemented mean and gaussian curvature methods in python and thought
> they may be appropriate for inclusion to skimage. An outline of the code is
> provided below. It assumes a surface defined as a function of two
> coordinates, e.g. z = Z(x, y). The curvature calculations are from the
> following two papers

Thank you for your post, and sorry for not responding earlier: I was
hoping someone more knowledgeable would jump in. Since I am not that
person, would you mind explaining to me the gist of this work, and
what it aims to do, perhaps with some example images?

Thanks!
Stéfan

[Brickle Macho]

Hi St�fan and others,

As requested.

Depth maps provide data on non-flat surfaces. Curvature is the amount
by which a object/surface deviates from being flat. Using differential
geometry, we can find the compute the curvature of the surfaces in an
image. Due to the structure of a depth map (essentially a Monge Patch
see: http://mathworld.wolfram.com/MongePatch.html), numerical
computation of the surface curvature map only requires the first and
second partial derivative estimates at each surface point and can be
estimated using np.gradient() function. Here is example output form the
provided code: http://i.imgur.com/kyMaR1F.jpg If interested, slightly
more detail provided after my signature block.

I have no ego over this, so if not appropriate to skimage that is ok.
Heck, someone else can take control and include this into skimage, but I
suspect like me, everyone is busy with life and other side projects.

Michael.
--


What follows is a more detailed version of above. Apologies if some of
the following is obvious some.

In my current research, my focus is on depth maps or range images.
Others in our research ground are using point clouds. The benefit of
depth maps is that they are basically a structured point cloud but allow
me to use common 2D image processing algorithms on the depth maps, and
not worry about a specific point cloud library or algorithms. For
example applying some simple morphological operations (from skimge) I
have been segment objects at varying depths, i.e segment objects that
standout.

A depth map is gray-scale image which encodes the distance for a certain
view or perspective. Depth maps record the first surface seen and
(obviously) cannot display information about occluded surfaces, surfaces
seen or refracted through transparent object or reflected in mirrors. I
also have an aligned with the color view of the same scene, but not
directly relevant to this post.

Curvature is the amount by which a object/surface deviates from being
flat. Curvature can be used to describe a surface. Either using the
curvature pixel value itself as the feature or calculating one of eight
surface primitives are: peak, pit, ridge, valley, flat, minimal, saddle
ridge, and saddle valley surfaces. In the code provided in the original
post I just calculate the the curvature value of each pixel, resulting
in a curvature feature map. Technically the code provided can be
applied to any grey-scale image, but really only makes sense when
applied to range iamges.

Depth maps provide data on non-flat surfaces. Differential geometry can
be used to process and analyse data on non-flat surfaces. Using
differential geometry, we can find the compute the curvature of the
surfaces in an image. A depth or range value at a point (x, y ) is
given by a single-valued function z(x,y), essentially a Monge Patch see:
http://mathworld.wolfram.com/MongePatch.html. Because of this the
Gaussian and mean curvatures are only related to the first and second
partial derivatives of the function z with respect to x and y. Thus
the numerical computation of the surface curvature map only requires the
first and second partial derivative estimates at each surface point.
These can be estimated using np.gradient() function. Once we have these
estimates, it is simply a mater of substituting into the equations from
the various papers.

Here is the code applied to a range image from the Middlebury data set
(http://vision.middlebury.edu/stereo/data/). The linked image shows he
original range image, mean and guassian curvature, and simplified mean
and guassian curvature. The RGB imagery is included for completeness but
is not used in the calculation. Curvature Example:
http://i.imgur.com/kyMaR1F.jpg

> St�fan
>