Average Curvature of a Closed 3D Surface

Hi, All --

I would like to compute the average curvature of a closed 3D surface for which I know the inside-outside function, F(x,y,z)=0. (FWIW, centered on the origin and symmetric about all three Cartesian axes.)

I have implemented a couple of different approaches for Gaussian and mean curvature using divergence, the Hessian, the first and second fundamental forms, and other approaches of which I have a similarly weak grasp. I am generally able to get things to work for degenrate cases (e.g., a sphere), but if I increase the complexity of my object, I am tending to get answers that imply I am flirting with the precipice of numerical instability (e.g., imaginary results where Re(x)/Im(x)~10^8).

I like the idea of this approach because it "feels clean": (1) use continuous function; (2) apply magic; (3) integrate in spherical coordinates; and (4) profit. However, I am beginning to think that I need to tune out the sires sweetly singing and go for a more practical (i.e., tractable) approach.

My other idea is to start by using the inside-outside function to generate [X,Y,Z] similar to what is used by surf(). I can then use surfature() from the FEX to get the curvatures at every point on my [X,Y,Z] mesh. Finally, I could convert everything to spherical coordinates and integrate the curvature arrays to get average values.

My concern here is that reregularly-spaced [X,Y,Z] does not translate to regularly-spaced [r,theta,phi], so the calculation is somehow pathologically biased. I am also unsure about how to do the actual integration after converting my array to spherical coordinates.

Apologies for such a long text-dense question. I thought about trying to illustrate this better with a couple of code snippets, but the algebra is really messy and not all that informative.

Thanks, in advance.

[Mathieu NOE]

hello

My concern here is that regularly-spaced [X,Y,Z] does not translate to regularly-spaced [r,theta,phi], so the calculation is somehow pathologically biased. I am also unsure about how to do the actual integration after converting my array to spherical coordinates.

if your integration process requires regularly-spaced [r,theta,phi], why not interpolate [r,theta,phi] obtained avec conversion from cartesian to spherical coordinates ?

[Moe Szyslak]

Hi, Mathieu -- thank you for the response. I misspoke in my original post. It's not that I need "regularly-spaced" data, but rather, that I would prefer the data to be uniformly distributed. Regularly-spaced [x,y,z] will result in [r,theta,phi] that are preferentially oriented at the poles. I suppose that I could convert using cart2sph() and then use those points to generate points that would be evenly distributed on the surface of my object, but that seems less-than-ideal.

In the end, though, with this approach I am still stuck with the fact that I will have an NxN array of curvatures that I am attempting to average over a surface defined by an NxNx3 array of coordinates, and I am not exactly sure how to do that.

[Mathieu NOE]

hello again

do you have kind of (simpler) working code and data to work with ?

[Moe Szyslak]

Yes, for sure. The discretization code (etc.) is attached. I think everything is pretty well-commented. I get to the point where I have my surface in spherical coordinates and my array of curvature values, but then how do I get a reliable average over all orientations? Thanks again.

[Mathieu NOE]

hello again

what do you mean by "a reliable average over all orientations?"

I suspect you want to dicretize the total surfaces into smaller areas (pavement) and do the integration (mean value of K ?) on these areas ?

[Paul]

Hi moe,

Can you provide a reference to the definiton of curvature for a surfact to be used for this problem?

[Moe Szyslak]

Hi, Paul -- thanks for th interest. Right now, I am working with Gaussian curvature, but honestly, I am somewhat agnostic about it. Please let me explain.

I have a collection of these closed surfaces and I am interested in quantifying the difference between their shapes. There are a number of metrics that can be used, with sphericity being arguably the most common. In 2D, roundness is a powerful metric, but it is somewhat difficult to quantify and not well-defined in 3D -- but I expect it to be related to curvature! So, that's what I am trying to do -- quantify curvature and compare it across surfaces. I would like to use the definition of curfature that provides the greatest distinction between shapes, and then try to figure out why that is the case.

Thanks again -- i would be happy if you had any thoughts on this.

[Paul]

Given that you have the form F(x,y,z) = 0, does eqn (16) in the MathWorld reference apply?

[Moe Szyslak]

Maybe. Indeed, this was my initial (and still preferred) approach. I think the catch is that what I have is actually the inside-outside function for the surface -- i.e., if F(x,y,z)=1 the point is on the surface, F<0 inside, and F>0 outside. I stated this clumsily (and incorrectly) in my original post.

There are several closed-form approaches that would appear to work (e.g., eq. 16), but all require a good bit of manipulation (differentiation, algebra) to get there. Then I am left with something that must be numerical integrated over some portion of the surface to get an average. I am OK with all of this, but I have struggled to obtain solutions that are free of numerical artifacts (e.g., very, very small imaginary parts), asymptotically reasonable, and entirely consistent with known solutions for simple cases (spheres and scalene ellipsoids). Tracing the problem (math error on my part? ill-posedness? incorrect selection/implementation of numerical tools? etc.?) has proved troublesome.

[Paul]

Seems that equation (16) still appliels to F(x,y,z) = 1 because that's the same as G(x,y,z) = 0 where G(x,y,z) = F(x,y,z) - 1, and all of the partials of F are the same as those of G.

Why not post a relatively simple example of an F(x,y,z) in which you're interested, or at least interested for testing?

[Moe Szyslak]

Paul -- sure thing! Here is the formula for my surfaces:

Without loss of generality for this calculation I assume a surface to be centered on the origin and aligned with the Cartesian axes, allowing us to remove the absolute values, which simplifies things somewhat. From this expression I have attempted to compute curvature over the positive octant of the surface using the Hessian approach and divergence formula (both, e.g., as outlined here) and also using the first and second fundamental forms, similar to the method oitlined in the MathWorld page linked above.

All three procedures are effectively equivalent, of course, and also equivalently messy because of the non-integer shape exponents (i.e., not much cancels out when you start taking derivatives and combining terms). It is clearly possible (likely, even) that I have made a mistake in any of the calculations. When I test my equations on degenerate cases I get correct solutions, but only sometimes (e.g., curvature of the unit sphere will be 1 at (0,0,1) but undefined at (1,1,0)). Results are equally frustrating for a scalene ellipsoid.

Thanks, in advance, for any thoughts that you have.

[Paul]

If you're only getting correct solutions sometimes ....

Here's an approach to computing the Gaussian curvature based on an equation in the link cited in the comment above.

Define a sphere

x = sym('x',[1 3]);
syms r real
F = sum(x.*x) - r^2

F = 

Compute the Gaussian curvature

gradF = gradient(F,x).';
hessF = hessian (F,x);

K_sphere = gradF*adjoint(hessF)*gradF.'/(gradF*gradF.')^2;
K_sphere(x) = simplify(subs(K_sphere,x(3)^2,rhs(isolate(F==0,x(3)^2))))

K_sphere(x1, x2, x3) = 

Define an ellipsoid and compute the Gaussian curvature

syms a b c positive
F = sum((x.^2./[a b c].^2)) - 1

F = 

gradF = gradient(F,x).';
hessF = hessian (F,x);

K_ellip = gradF*adjoint(hessF)*gradF.'/(gradF*gradF.')^2;
K_ellip(x) = simplify(subs(K_ellip,x(3)^2,rhs(isolate(F==0,x(3)^2))))

K_ellip(x1, x2, x3) = 

The result matches eqn (16) at Ellipsoid

Both results above took an extra step to get K as a function of two or fewer variables.

Try with the equation above, ignoring absolute values as directed ...

syms epsilon_1 epsilon_2 positive
r = sym('r',[1 3],'positive');

F = sum((x(1:2)./r(1:2)).^(2/epsilon_2))^(epsilon_2/epsilon_1) + (x(3)/r(3))^(2/epsilon_1) - 1

F = 

gradF = gradient(F,x).';
hessF = hessian (F,x);

K_moe = gradF*adjoint(hessF)*gradF.'/(gradF*gradF.')^2;
K_moe = simplify(K_moe)

K_moe = 

If neccesary to eliminate one of the variables, the best I could do for now is (ignoring analytic constraints always is a bit uncomfortable because I don't know why that directive is needed)

isolate(F==0,(x(3)/r(3))^(2/epsilon_1))

ans = 

sol = solve(ans,x(3),'IgnoreAnalyticConstraints',true,'ReturnConditions',true)

sol = struct with fields:
            x3: r3*(1 - ((x1/r1)^(2/epsilon_2) + (x2/r2)^(2/epsilon_2))^(epsilon_2/epsilon_1))^(epsilon_1/2)
    parameters: [1x0 sym]
    conditions: symtrue

sol.x3

ans = 

K_moe(x) = subs(K_moe,x(3),sol.x3)

K_moe(x1, x2, x3) = 

That's not a very nice expression, but it should be correct unless the original equation for the curvature doesn't apply to this particular surface or if that solve for x(3) isn't true for some reason.

Here's one sanity check that seems correct

subs(K_moe,[epsilon_1,epsilon_2],[2,2])

ans(x1, x2, x3) = 

0

At this point, I'm not sure what you'd want to do with any of the K expressions. I gather you want to compute some sort of "average" curvature as a figure of merit. How would that be defined? As a surface integral of K over the surface normalized by something?

[Moe Szyslak]

Hey, Paul -- this is fantastic -- thanks so much. This is exactly what I was trying to do in developing a closed-form solution to the problem.

I am not very familiar with MATLAB's symbolic math functionality. and I think where I may have been going wrong was in attempting to develop the analytical solution independently (using Mathcad, gasp!) prior to implementing in MATLAB. Your approach is much cleaner and more robust. I am going to play around with this some and see if I can either break it or convince myself that it does what I think it does.

Ultimately, yes, I am looking to get the average value of curvature over the entire surface. My instinct is to apply a spherical transformation and then integrate over the range (exploiting symmetry). My concern with this, however, is that all of the "action" is happening in areas that are near one of the Cartesian coordinate planes -- that is, at the bounds of my integration. I am wondering if it makes more sense to take the average of the absolute value of the curvature over the entire surface () rather than just 1/8th of it.

Thanks again for all of your help. I am learning a lot.

Answer 1

[Mathieu NOE]

that's me again ! :)

so a bit for my own fun , I tried a few things

here I decided I wanted to define the intersection of your surface with a 3D rectangular box (any other shape can be done , and maybe here the spherical coordinates can be useful to create other pavement shapes)

from the "red dots" carpet you can compute the integral (or mean) of K (or H) over this area with the triangulation method

you can reduce the constrain of using a very refined mesh as the triangulation method is pretty accurate even with fairly coarse resolution.

if you compare the result of the triangulation mean with the simple average values of selected K values (w/o using the uneven areas info) you see both methods will not give you the same results (that's obvious , sorry ! ) , unless you use a very refined mesh resolution, but that is not my recommendation.

in other words , keep the resolution at reasonnable value (N = 100 seems to suffice in my eyes) and use the triangulation integration method (the rest is just fyi)

hope it helps !

% Driver script to create a superquadric, calculate the surfavece area and
% volume, generate a surface mesh in Cartesian coordinates, compute the
% curvature, and convert the surface mesh to spherical coordinates.

% Inside-outside function for superquadric
% f(x,r,eps) =
% ((x/rx)^(2/eps2)+(y/ry)^(2/eps2))^(eps2/eps1)+(z/rz)^(2/eps1)


% Preliminaries
close all;
clear;

% Input Values
rx = 0.6; ry = 1.2; rz = 1.0;       % semiaxis lengths
eps1 = 1.4; eps2 = 0.6;             % shape exponents
n = 100;                           % resolution for discretization

% Compute volume and surface area
V = 2*rx.*ry.*rz.*eps1.*eps2.*beta(eps1/2,eps1+1).*beta(eps2/2,(eps2+2)/2);
SA = superellipsoid_SA2(rx,ry,rz,eps1,eps2);

% Discretize and compute curvature
[X,Y,Z] = superquad(rx,ry,rz,eps1,eps2,n);
[K,H,P1,P2] = surfature(X,Y,Z);

% Plot object and colormap of curvature
surfl(X,Y,Z), shading interp, colormap(copper), axis equal
figure;
% surf(X,Y,Z,K,'facecolor','interp','edgecolor','none');
surf(X,Y,Z,K);
set(gca,'clim',[-1,1]);
colorbar;
axis equal;

% % Convert to spherical coordinates
% [az,el,rho] = cart2sph(X,Y,Z);

% extract a local area (intersection of you surface with a rectangular box)
indX = X>-0.5 & X < -0.3;
indY = Y>-0.75 & Y < 0;
indZ = Z> 0.2 & Z < 0.5;
ind = indX & indY & indZ;
xlocal = X(ind);
ylocal = Y(ind);
zlocal = Z(ind);

hold on 
scatter3(xlocal,ylocal,zlocal,25,"r","filled");

%% 2D integration using triangulation:
% Let's say your integration area is given by the triangulation "trep", 
% which for example could be obtained by trep = delaunayTriangulation(x(:), y(:)).
% If you have your values z corresponding to z(i) = f(trep.Points(i,1), trep.Points(i,2)), 
% you can use the following integration routine. 
% It computes the exact integral of the linear interpolant. 
% This is done by evaluating the areas of all the triangles and then using these areas
% as weights for the midpoint(mean)-value on each triangle.

trep = delaunayTriangulation(xlocal(:), ylocal(:));

%% THIS code below can be deleted latter on - just fyi so you don't have to ask you the same question as I did
% xtrep = trep.Points(:,1);
% ytrep = trep.Points(:,2);
% % NB xtrep = xlocal and same for ytrep = ylocal. No need to create interpolated z values here
% % proof in the plot below
% figure,
% plot(xlocal,ylocal,'*r');
% hold on
% plot(xtrep,ytrep,'ob');

% let's do the integration now 
[intZ,meanZ] = integrateTriangulation(trep, zlocal)

% now compare meanZ (true mean) with mean(zlocal) (is of course a crude approx as we don't take
% into account the non uniform grid into consideration)
zlocal_mean = mean(zlocal)
% by refining the mesh , the two values will converge . 
% The triangulation method allows a better and faster approach
% as you can work with a coarser mesh without loss of accuracy

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [intZ,meanZ] = integrateTriangulation(trep, z)
    P = trep.Points; T = trep.ConnectivityList;
    d21 = P(T(:,2),:)-P(T(:,1),:);
    d31 = P(T(:,3),:)-P(T(:,1),:);
    areas = abs(1/2*(d21(:,1).*d31(:,2)-d21(:,2).*d31(:,1)));
    % integral of Z
    intZ = areas'*mean(z(T),2);
    % mean of Z = intZ divided by total area
    meanZ = intZ/(trapz(areas));
end

[Moe Szyslak]

Thank you so much for the continued help. I think that I mostly understand what you've done, but I am not certain that it gets me where I am trying to go, so I have a couple of more questions.

Since I am trying to get the average curvature over the whole surface and the object is 3D symmetric, I am using the following bounds:

indX = X >= 0 & X <= max(X);

indY = Y >= 0 & Y <= max(Y);

indZ = Z >= 0 & Z <= max(Z);

This defines a curved surface, but I think that this command produces a planar triangulation:

trep = delaunayTriangulation(xlocal(:), ylocal(:));

Then this function computes the average z-coordinate for each triangle and, subsequently, the area-weighted average z-coordinate for the whole surface:

[intZ,meanZ] = integrateTriangulation(trep, zlocal);

Is my understanding correct so far? Assuming so, I have a couple of additional questions.

Because the triangulation is planar, the areas of the triangles are not the same as the "curvilinear" triangles that would result if the triangulation were carried out on the surface. So, I tried to do the triangulation in 3D, but as I am sure you already know, delaunayTriangulation(xlocal,ylocal,zlocal) returns the volumetric triangulation, not the surface triangulation that I desire. Is there a way to just get the surface triangulation?

Once I have the 3D surface triangulation, I assume that I could use a slightly modified version of your approach to get the average curvature for each triangular patch, which could then be used to compute the average curvature for the whole surface. Correct?

A final question -- why do you use trapz(areas) to get the total area instead of sum(areas)? The results are obviously very close, but not the same.

Thanks again.

[Moe Szyslak]

OK, thanks to your help, Mathieu, I think that the script below does mostly what I want.

My concern is that I observe asymptotic behavior only for very high resolutions. For reference, my laptop (16 GB RAM) protests n=10000 and does n=5000 in about 4 min. My workstation (64 GB RAM) takes about 5 min for n=10000 and 15 min for n = 15000, but throws an error for n=20000. So, we're getting to the outer limits of what is reasonable, given that I have many thousands of surfaces to characterize.

n meanK

100 1.203

200 1.295

500 1.423

1000 1.528

2000 1.642

5000 1.817

10000 1.974

15000 2.076

% Driver script to create a superquadric, calculate the surfavece area and
% volume, generate a surface mesh in Cartesian coordinates, compute the
% curvature, and convert the surface mesh to spherical coordinates.
%
% Acknowledgements:
% Significant assistance from Mathieu Noe of the MATLAB community.
% www.mathworks.com/matlabcentral/answers/2172786-average-curvature-of-a-closed-3d-surface

% Inside-outside function for superquadric
% f(x,r,eps) =
% ((x/rx)^(2/eps2)+(y/ry)^(2/eps2))^(eps2/eps1)+(z/rz)^(2/eps1)

% Preliminaries
close all;
clear;

% Input Values
rx = 1.25; ry = 0.75; rz = 1.0;     % semiaxis lengths
eps1 = 1.5; eps2 = 0.5;             % shape exponents
n = 1000;                           % resolution for discretization

% Compute volume and surface area
V = 2*rx.*ry.*rz.*eps1.*eps2.*beta(eps1/2,eps1+1).*beta(eps2/2,(eps2+2)/2);
SA = superellipsoid_SA2(rx,ry,rz,eps1,eps2);

% Discretize and compute curvature
[X,Y,Z] = superquad(rx,ry,rz,eps1,eps2,n);
[K,H,P1,P2] = surfature(X,Y,Z);

% Plot object and colormap of curvature
%surfl(X,Y,Z), shading interp, colormap(copper), axis equal
figure;
surf(X,Y,Z,K,'facecolor','interp','edgecolor','none');
set(gca,'clim',[-5,5]);
colorbar;
axis equal;

% extract a local area (top half of the surface)
indX = X >= min(min(X)) & X <= max(max(X));
indY = Y >= min(min(Y)) & Y <= max(max(Y));
indZ = Z >= 0           & Z <= max(max(Z));
indK = ~isnan(K);
ind = indX & indY & indZ & indK;

A = unique([X(ind) Y(ind) Z(ind) K(ind)],'rows','stable');

xlocal = A(:,1);
ylocal = A(:,2);
zlocal = A(:,3);
Klocal = A(:,4);

hold on 
scatter3(xlocal,ylocal,zlocal,25,"r","filled");

% 2D integration using triangulation:
trep = delaunayTriangulation(xlocal(:), ylocal(:));
[intK,meanK,sArea] = integrateTriangulation(trep, zlocal, Klocal);

% COmpute naive average curvature
Klocal_mean = mean(Klocal);

% Output results
fprintf('Surface area (closed form):\t\t%f\n', SA);
fprintf('Surface area (triangulation):\t%f\n', 2*sArea);
fprintf('Area-weighted curvature:\t\t%f\n', meanK);
fprintf('Naive average curvature:\t\t%f\n', Klocal_mean);
fprintf('Particle volume:\t\t\t\t%f\n', V);

%% Support Functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [intK,meanK,sArea] = integrateTriangulation(trep, z, K)
    P = trep.Points; T = trep.ConnectivityList;
    V = [P z];
    d21 = V(T(:,2),:)-V(T(:,1),:);
    d31 = V(T(:,3),:)-V(T(:,1),:);
    areas = vecnorm(cross(d21,d31)')/2;
    % areas = abs(1/2*(d21(:,1).*d31(:,2)-d21(:,2).*d31(:,1)));
    % integral of Z
    intK = areas*mean(K(T),2);
    % mean of Z = intZ divided by total area
    sArea = sum(areas);
    meanK = intK/sArea;
end

%eof

Anybody have any thoughts about how to get a "true" solution without so much resolution dependency?

[Mathieu NOE]

hello again (and tx for accepting my suggestions ! ) :)

I am not sure that is needed indeed :

indK = ~isnan(K);

ind = indX & indY & indZ & indK;

A = unique([X(ind) Y(ind) Z(ind) K(ind)],'rows','stable');

xlocal = A(:,1);

ylocal = A(:,2);

zlocal = A(:,3);

Klocal = A(:,4);

but what I see as possible issue in the existing code is that surfature works for regularly gridded surfaces with is not the case of your superquadric (see the way it is parametrized inside the superquad function)

so I made a test between the regular code (SQ_curvDriver_loop.m) and a modified one with linear interpolation using scatteredInterpolant (SQ_curvDriver_loop1.m)

using the same boundary conditions as you , but slightly corrected (in "modern" style) :

indX = X >= min(X,[],'all') & X <= max(X,[],'all');
indY = Y >= min(Y,[],'all') & Y <= max(Y,[],'all');
indZ = Z >= 0           & Z <= max(Z,[],'all');

the regular code (without any other modification than running a for loop for n = 100 300 700 1500 3100) will generate this result :

the modified code with scatteredInterpolant will generate this result :

seems that now we have a feeling that both methods converge , and to the same limit !!

what would be great is now to have the "true" analytical result to compare with...

now if you want to avoid the additionnal burden of scatteredInterpolant , you need to find a way to change the way the surface is parametrized and use instead linear gridded x, y data.

both codes are in attachment

hope it helps !