Thread Subject:

Curvature of a 3D implicit function

Hi,

Given an implicit function f (x,y,z) = 0 , is there a formula that I can use to calculate the curvature at a point?

Yours,
Weng

"Weng Boon" <chiaweng@gmail.com> wrote in message <ics6dc$khv$1@fred.mathworks.com>...
> Hi,
>
> Given an implicit function f (x,y,z) = 0 , is there a formula that I can use to calculate the curvature at a point?
>
> Yours,
> Weng
- - - - - - - - - -
  Your implicit function will define a surface in three dimensions. Curvature at a point on a surface is only defined in terms of curves that follow the surface through the point. The maximum and the minimum such curvatures on a surface are called its "principal curvatures". There are some methods of determining these for a point on a surface in terms of various of its derivatives.

  You might get a start on the subject by looking at the Wikepedia website at:

 http://en.wikipedia.org/wiki/Curvature

It has a number of links to further discussions of the subject. This subject is taught in differential geometry classes at universities.

  I'm afraid my last class in differential geometry was some six decades ago, so the subject would require some extensive work on my part before I could discuss it intelligently.

Roger Stafford

"Weng Boon" <chiaweng@gmail.com> wrote in message <ics6dc$khv$1@fred.mathworks.com>...
> Hi,
>
> Given an implicit function f (x,y,z) = 0 , is there a formula that I can use to calculate the curvature at a point?

At any point on the surface, The implicit can be easily transformed to explicit form

Z = F(X,Y)

simply by choosing the orthonormal basis such that the Z coordinate parallel to gradient f at the point under consideration.

Then take the Hessian of F, the two eigen values of the Hessian are the two principal curvatures.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ict53u$bdh$1@fred.mathworks.com>...
> "Weng Boon" <chiaweng@gmail.com> wrote in message <ics6dc$khv$1@fred.mathworks.com>...
> > Hi,
> >
> > Given an implicit function f (x,y,z) = 0 , is there a formula that I can use to calculate the curvature at a point?
>
> At any point on the surface, The implicit can be easily transformed to explicit form
>
> Z = F(X,Y)
>
> simply by choosing the orthonormal basis such that the Z coordinate parallel to gradient f at the point under consideration.
>
> Then take the Hessian of F, the two eigen values of the Hessian are the two principal curvatures.
>
> Bruno
- - - - - - - - -
  Bruno has given you a good answer, Weng. I would like to propose a modification of his method however which would not require conversion from the implicit equation f(x,y,z)=0 to an explicit form such as his Z=F(X,Y) in case that would be more convenient for you.

  If the point (x0,y0,z0) on the surface is the one at which the principal curvatures are to be determined, let fx, fy, fz, fxx, fxy, fxz, fyx, fyy, fyz, fzx, fzy, and fzz, be the three first and the nine second partial derivatives of f evaluated at this point.

  In the absence of an explicit form it is necessary to decide which side of the surface defined by f(x,y,z) = 0 is to be regarded as the side the surface bends to in order to be regarded as positive curvature. If f is positive on the side of positive curvature, the 's' quantity below should be set at +1. Otherwise it should be set to -1.

 s = +1; % Positive curvature on positive side of f, otherwise set s = -1
 N = [fx;fy;fz]; % The gradient of f at (x0,y0,z0)
 M = null(N.'); % Vectors in M and the gradient N are mutually orthogonal
 H = [fxx,fxy,fxz;fyx,fyy,fyz;fzx,fzy,fzz]; % The Hessian of f at (x0,y0,z0)
 [U,D] = eig(M.'*H*M); % These eigenvalues are proportional to curvatures
 K = -s*diag(D).'/norm(N); % Convert to the true curvatures
 V = M*U; % Obtain the corresponding eigenvectors in curvature directions

  The two-element row vector K will have the two principal curvatures, and the two corresponding columns of V will each be unit orthogonal vectors in the two directions along which those principal curvatures occur.

  I should point out that these curvatures are those given by curves on the surface through the given point which lie in planes containing the normal to the surface at that point. Other such planar curves through the point in intermediate directions will have curvatures between these two extreme values.
 
Roger Stafford

Good work out by Roger to derive the formula that can be directly applied to implicit surface.

Bruno

On Dec 1, 8:13 pm, "Roger Stafford"
<ellieandrogerxy...@mindspring.com.invalid> wrote:
>   In the absence of an explicit form it is necessary to decide which side of the surface defined by f(x,y,z) = 0 is to be regarded as the side the surface bends to in order to be regarded as positive curvature.  If f is positive on the side of positive curvature, the 's' quantity below should be set at +1.  Otherwise it should be set to -1.
>
>  s = +1; % Positive curvature on positive side of f, otherwise set s = -1
>  N = [fx;fy;fz]; % The gradient of f at (x0,y0,z0)
>  M = null(N.'); % Vectors in M and the gradient N are mutually orthogonal
>  H = [fxx,fxy,fxz;fyx,fyy,fyz;fzx,fzy,fzz]; % The Hessian of f at (x0,y0,z0)
>  [U,D] = eig(M.'*H*M); % These eigenvalues are proportional to curvatures
>  K = -s*diag(D).'/norm(N); % Convert to the true curvatures
>  V = M*U; % Obtain the corresponding eigenvectors in curvature directions
>
>   The two-element row vector K will have the two principal curvatures, and the two corresponding columns of V will each be unit orthogonal vectors in the two directions along which those principal curvatures occur.
>
>   I should point out that these curvatures are those given by curves on the surface through the given point which lie in planes containing the normal to the surface at that point.  Other such planar curves through the point in intermediate directions will have curvatures between these two extreme values.
>



Thanks so much, Roger,
Could you recommend a textbook or reference where this derivation is
discussed in some more details?
Igor

rych <rychphd@gmail.com> wrote in message <6ed0b534-62a3-4cff-bd7a-c0990bdd500e@h17g2000pre.googlegroups.com>...
> ......
> Thanks so much, Roger,
> Could you recommend a textbook or reference where this derivation is
> discussed in some more details?
> Igor
- - - - - - - - -
  I am not good at giving references to books or papers on this subject, since my instruction in that area occurred so long ago, but I can make an attempt to justify the algorithm I presented.

  Starting with an implicitly defined surface, f(x,y,z) = 0, we wish to find the curvatures at some point (x0,y0,z0) on this surface. As Bruno asserts, in principle we should be able to find a rotated orthogonal coordinate system (u,v,w) such that w = h(u,v) would identically satisfy

 f(u,v,h(u,v)) = 0

in the neighborhood of (x0,y0,z0) and such that at the point itself the gradient of f should be parallel to the w-axis. The orientation of u and v around the w-axis would be arbitrary.

  We therefore define

 F(u,v) = f(u,v,h(u,v))

for some such presumed system. Since F is identically zero in the neighborhood, all its derivatives are also zero, so we have

 Fu = fu + fw*hu = 0
 Fv = fv + fw*hv = 0

where Fu, fu, and hu indicate partial derivatives with respect to u and similarly for v and w. Also we have

 Fuu = (fuu + fuw*hu) + (fwu + fww*hu)*hu + fw*huu = 0
 Fuv = (fuv + fuw*hv) + (fwv + fww*hv)*hu + fw*huv = 0
 Fvu = (fvu + fvw*hu) + (fwu + fww*hu)*hv + fw*hvu = 0
 Fvv = (fvv + fvw*hv) + (fwv + fww*hv)*hv + fw*hvv = 0

At the point (x0,y0,z0) by the above assumption we must have fu = fy = 0, which means that hu and hv are also zero at the point and thus that

 fuu + fw*huu = 0
 fuv + fw*huv = 0
 fvu + fw*hvu = 0
 fvv + fw*hvv = 0

at the point. This means that the eigenvalues of the 2 x 2 Hessian matrix

 [huu huv
  hvu hvv]

which give the principal curvatures will be equal to -1/fw times those of the Hessian matrix

 [fuu fuv
  fvu fvv]

at the point. By this means we have thus avoided specifically finding w = h(u,v). Finding the curvatures does not require it.

  In the code I presented, the magnitude of the gradient N (norm(N)) is this quantity fw, provided the w-axis is oriented in the same direction as the gradient. The product M.'*H*M converts the 3 x 3 Hessian matrix of f(x,y,z) to the 2 x 2 Hessian of [fuu,fuv;fvu,fvv] for the particular orientation of the u and v axes selected by taking the 'null' of N.

  In case w is oriented opposite to the gradient of f, this means that the curvature is regarded as being positive if the curvature is in the opposite direction of the gradient and hence the 's' quantity would be -1 to compensate.

  As for V it is necessary to rotate the eigenvectors relative to the u,v axes back to those of x, y, and z axes in the expression V = M*U. These adjusted eigenvectors will then give the directions relative to the x, y, z axes of the principal curvatures.

  I neglected to mention earlier that curvatures along planar curves in intermediate directions are given by K1*cos(theta)^2+K2*sin(theta)^2 where K1 and K2 are the respective curvatures and theta the angle between the intermediate direction and the K1 direction.

  I hope all this at least gives you a rough idea of the reasoning involved.

Roger Stafford

Thanks Bruno and Roger! I apologise for my slow response. I was on a holiday. Roger, it works! I appreciate it.

Weng