Dear all,
Please advise in which module can I find implementation of Gaussian Curvature of a surface?
Something similar to the MATLAB function: Surface Curvature - File Exchange - MATLAB Central
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.
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!
My new library Comodo can do curvature analysis on arbitrary surfaces, see this demo:
1 Like