# Gaussian Curvature calculation

Dear all,

Please advise in which module can I find implementation of Gaussian Curvature of a surface?

Something similar to the MATLAB function: https://www.mathworks.com/matlabcentral/fileexchange/11168-surface-curvature

Welcome tho the Julialang Discourse!

For what i saw, there isn’t a package to do exactly what you want, with that said, That routine doesn’t seem too difficult to port to Julia

1 Like

Ok…being newbie for Julia I have to admit that your DIY advise was the most correct way of doing things )

My implementation is below. `gradient` and `dot` utilities are customized to my 2d case.
I totally agree it’s not generalize and probably looks ugly for experienced Julia devs but…it works )

``````using LinearAlgebra

"Numerical `gradient` calculation over 1st matrix axis Y (rows)"
dv  = diff(v; dims=1)/2
a   = [dv[,:]; dv]
a .+= [dv; dv[[end],:]]
a
end

"Numerical `gradient` calculation over 2nd matrix axis X (cols)"
dv  = diff(v; dims=2)/2
a   = [dv[:,] dv]
a .+= [dv dv[:,[end]]]
a
end

"""
Numerical gradient of a matrix over first 2 axes - Y and X

Output corresponds to matlab `gradient` function.
First output parameter is gradient matrix for 2nd axis - X. Second output is for 1st axis - Y
"""
end

"""
Calculate dot product over matrices 1st axis Y (rows)

Returns: vector with elements which are dot product of corresponding matrices rows
"""
function doty(a::AbstractMatrix, b::AbstractMatrix)
[dot(a[i,:], b[i,:]) for i in 1:size(a)]
end

"""
Calculate dot product over matrices 2nd axis X (cols)

Returns: vector with elements which are dot product of corresponding matrices rows
"""
function dotx(a::AbstractMatrix, b::AbstractMatrix)
[dot(a[:,i], b[:,i]) for i in 1:size(a)]'
end

"""
Calculate cross product over matrices 1st axis Y (rows)

Returns: vector with elements which are cross product of corresponding matrices rows
"""
function crossy(a::AbstractMatrix, b::AbstractMatrix)
c = [cross(a[i,:], b[i,:]) for i in 1:size(a)]
Matrix(hcat(c...)')
end

"""
Calculate cross product over matrices 2nd axis X (cols)

Returns: vector with elements which are cross product of corresponding matrices rows
"""
function crossx(a::AbstractMatrix, b::AbstractMatrix)
c = [cross(a[:,i], b[:,i]) for i in 1:size(a)]'
Matrix(vcat(c...)')
end

"""
Gaussian, mean, min and max curvatures of a surface

Gaussian and Mean curvatures
`k,h = surfature(x,y,z)`, where x,y,z are 2d arrays of points on the surface.
k and h are the gaussian and mean curvatures, respectively.

`surfature` returns 2 additional arguments: `k,h,pmax,pmin = surfature(x,y,z)`.
pmax and pmin are the minimum and maximum curvatures at each point, respectively.

Function is shamelessly plagiated from matlab one.
See https://www.mathworks.com/matlabcentral/fileexchange/11168-surface-curvature
"""
function surfature(X::AbstractMatrix, Y::AbstractMatrix, Z::AbstractMatrix)

# First Derivatives

# Second Derivatives

# Reshape 2D Arrays into Vectors
Xu = Xu[:];   Yu = Yu[:];   Zu = Zu[:];
Xv = Xv[:];   Yv = Yv[:];   Zv = Zv[:];
Xuu = Xuu[:]; Yuu = Yuu[:]; Zuu = Zuu[:];
Xuv = Xuv[:]; Yuv = Yuv[:]; Zuv = Zuv[:];
Xvv = Xvv[:]; Yvv = Yvv[:]; Zvv = Zvv[:];

Xu          =   [Xu Yu Zu]
Xv          =   [Xv Yv Zv]
Xuu         =   [Xuu Yuu Zuu]
Xuv         =   [Xuv Yuv Zuv]
Xvv         =   [Xvv Yvv Zvv]

# First fundamental Coeffecients of the surface (E,F,G)
E           =   doty(Xu,Xu)
F           =   doty(Xu,Xv)
G           =   doty(Xv,Xv)

m           =   crossy(Xu,Xv)
p           =   sqrt.(doty(m,m))
n           =   m./[p p p]

# Second fundamental Coeffecients of the surface (L,M,N)
L           =   doty(Xuu,n)
M           =   doty(Xuv,n)
N           =   doty(Xvv,n)

s, t = size(Z)

# Gaussian Curvature
K = (L .* N - M .^ 2) ./ (E .* G - F .^ 2)
K = reshape(K,s,t)

# Mean Curvature
H = (E .* N + G .* L - 2 .* F .* M) ./ (2 * (E .* G - F .^ 2))
H = reshape(H,s,t)

# Principal Curvatures
Pmax = H + sqrt.(H .^ 2 - K)
Pmin = H - sqrt.(H .^ 2 - K)

K, H, Pmax, Pmin
end
``````
1 Like

Do you mean a surface represented by an equation z = f(x, y)?

Please do not post code that does not have an open source license.

@dpsanders

Hm…but it also doesn’t have any license at all as I see

And it was posted on official exchange matlab service.

Here its python implementation posted on SO

So, I’m thinking it’s free for use

Yes, function of 2 variables `z = f(x,y)`

I mean do you have the function f explicitly?
If so, you can use ForwardDiff.jl to calculate the derivatives.

Ah…no, no explicit form. Just 2d arrays from image files
Thanks for pointing out `ForwardDiff.jl`

1 Like

@Ameba_Brain
Hi there, I am trying to calculate the different curvatures at each point in point-cloud. I have yet to find a library or an implementation in Julia other than yours.
Did you implement this code on point clouds by any chance? if so, can you please share an example?