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

Thank you in advance

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)"
function gradienty(v::AbstractMatrix)
    dv  = diff(v; dims=1)/2
    a   = [dv[[1],:]; dv]
    a .+= [dv; dv[[end],:]]
    a
end

"Numerical `gradient` calculation over 2nd matrix axis X (cols)"
function gradientx(v::AbstractMatrix)
    dv  = diff(v; dims=2)/2
    a   = [dv[:,[1]] 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
"""
function gradient(v::AbstractMatrix)
    return gradientx(v), gradienty(v)
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)[1]]
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)[2]]'
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)[1]]
    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)[2]]'
    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
    Xu,Xv = gradient(X)
    Yu,Yv = gradient(Y)
    Zu,Zv = gradient(Z)

    # Second Derivatives
    Xuu,Xuv = gradient(Xu)
    Yuu,Yuv = gradient(Yu)
    Zuu,Zuv = gradient(Zu)

    Xuv,Xvv = gradient(Xv)
    Yuv,Yvv = gradient(Yv)
    Zuv,Zvv = gradient(Zv)

    # 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?
Thanks in advance!