From 053a822221c8484b8fd1361ee177159f6ea5bebc Mon Sep 17 00:00:00 2001 From: Kevin-Mattheus-Moerman Date: Mon, 1 Apr 2024 22:49:20 +0100 Subject: [PATCH] Added tests for mesh_curvature_polynomial #20 --- examples/demo_mesh_curvature_polynomial.jl | 41 +++++++++------ src/functions.jl | 2 +- test/runtests.jl | 59 ++++++++++++++++++++++ 3 files changed, 86 insertions(+), 16 deletions(-) diff --git a/examples/demo_mesh_curvature_polynomial.jl b/examples/demo_mesh_curvature_polynomial.jl index 55f81e1..ae0877c 100644 --- a/examples/demo_mesh_curvature_polynomial.jl +++ b/examples/demo_mesh_curvature_polynomial.jl @@ -4,21 +4,30 @@ using GeometryBasics using FileIO # Example geometry -testCase = 6 +testCase = 7 if testCase == 1 r = 25.25 F,V = geosphere(5,r) # V = [GeometryBasics.Point{3, Float64}(v[1],v[2],3*v[3]) for v ∈ V] elseif testCase==2 - r = 1.0 + r = 25.25 F,V = quadsphere(3,r) elseif testCase==3 - r=sqrt(3) + r = sqrt(3) M = cube(r) - F=faces(M) - V=coordinates(M) + F = faces(M) + V = coordinates(M) # F = quad2tri(F,V; convert_method = "angle") -elseif testCase==4 # Merged STL for single object +elseif testCase==4 + r = 25 + s = r/10 + nc = round(Int64,(2*pi*r)/s) + d = r*2 + Vc = circlepoints(r, nc; dir=:cw) + num_steps = round(Int64,d/s) + num_steps = num_steps + Int64(iseven(num_steps)) + F, V = extrudecurve(Vc, d; s=1, n=[0.0,0.0,1.0], num_steps=num_steps, close_loop=true, face_type=:quad) +elseif testCase==5 # Merged STL for single object # Loading a mesh fileName_mesh = joinpath(comododir(),"assets","stl","stanford_bunny_low.stl") M = load(fileName_mesh) @@ -28,7 +37,7 @@ elseif testCase==4 # Merged STL for single object V = togeometrybasics_points(coordinates(M)) F,V,_ = mergevertices(F,V) F,V = subtri(F,V,2; method = :loop) -elseif testCase==5 # Merged STL for single object +elseif testCase==6 # Merged STL for single object # Loading a mesh fileName_mesh = joinpath(comododir(),"assets","stl","david.stl") M = load(fileName_mesh) @@ -37,7 +46,7 @@ elseif testCase==5 # Merged STL for single object F = togeometrybasics_faces(faces(M)) V = togeometrybasics_points(coordinates(M)) F,V,_ = mergevertices(F,V) -elseif testCase==6 +elseif testCase==7 # Loading a mesh fileName_mesh = joinpath(comododir(),"assets","obj","motherChild_5k.obj") M = load(fileName_mesh) @@ -52,26 +61,28 @@ M = GeometryBasics.Mesh(V,F) K1,K2,U1,U2,H,G = mesh_curvature_polynomial(F,V) ## Visualization - +s = pointspacingmean(F,V) cMap = :Spectral fig = Figure(size=(800,800)) ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "1st principal curvature") -hp1 =mesh!(ax1,M, color=K1, shading = FastShading, transparency=false,colormap=cMap,colorrange = maximum(abs.(K1)).*0.1.*(-1,1)) -hpn1 = dirplot(ax1,V,U1; color=:black,linewidth=1,scaleval=3,style=:through) +hp1 = mesh!(ax1,M, color=K1, shading = FastShading, transparency=false,colormap=cMap,colorrange = maximum(abs.(K1)).*0.1.*(-1,1)) +# hp1 = poly!(ax1,M, color=K1, shading = FastShading, transparency=false,colormap=cMap,colorrange = maximum(abs.(K1)).*0.1.*(-1,1),strokecolor=:black,strokewidth=2) +# normalplot(ax1,M) +hpn1 = dirplot(ax1,V,U1; color=:black,linewidth=2,scaleval=s,style=:through) Colorbar(fig[1, 2],hp1) ax1 = Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "2nd principal curvature") -hp2 =mesh!(ax1,M, color=K2, shading = FastShading, transparency=false,colormap=cMap,colorrange = maximum(abs.(K2)).*0.1.*(-1,1)) -hpn2 = dirplot(ax1,V,U2; color=:black,linewidth=1,scaleval=3,style=:through) +hp2 = mesh!(ax1,M, color=K2, shading = FastShading, transparency=false,colormap=cMap,colorrange = maximum(abs.(K2)).*0.1.*(-1,1)) +hpn2 = dirplot(ax1,V,U2; color=:black,linewidth=2,scaleval=s,style=:through) Colorbar(fig[1, 4],hp1) ax1 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Mean curvature") -hp2 =mesh!(ax1,M, color=H, shading = FastShading, transparency=false,colormap=cMap,colorrange = maximum(abs.(H)).*0.1.*(-1,1)) +hp2 = mesh!(ax1,M, color=H, shading = FastShading, transparency=false,colormap=cMap,colorrange = maximum(abs.(H)).*0.1.*(-1,1)) Colorbar(fig[2, 2],hp1) ax1 = Axis3(fig[2, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Gaussian curvature") -hp2 =mesh!(ax1,M, color=G, shading = FastShading, transparency=false,colormap=cMap,colorrange = maximum(abs.(G)).*0.025.*(-1,1)) +hp2 = mesh!(ax1,M, color=G, shading = FastShading, transparency=false,colormap=cMap,colorrange = maximum(abs.(G)).*0.025.*(-1,1)) Colorbar(fig[2, 4],hp1) fig diff --git a/src/functions.jl b/src/functions.jl index c79caa4..0126d83 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -2445,7 +2445,7 @@ which features a nice overview of the theory/steps involved in this algorithm. # References [F. Cazals and M. Pouget, "Estimating differential quantities using polynomial fitting of osculating jets", Computer Aided Geometric Design, vol. 22, no. 2, pp. 121-146, Feb. 2005, doi: 10.1016/j.cagd.2004.09.004](https://doi.org/10.1016/j.cagd.2004.09.004) """ -function mesh_curvature_polynomial(F::Vector{TriangleFace{Int64}},V::Vector{Point3{Float64}}) +function mesh_curvature_polynomial(F::Array{NgonFace{M, Int64}, 1},V::Vector{Point3{Float64}}) where M # Get the unique mesh edges E_uni = meshedges(F;unique_only=true) diff --git a/test/runtests.jl b/test/runtests.jl index 7e58883..8e05716 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2129,6 +2129,65 @@ end end +@testset "mesh_curvature_polynomial" verbose = true begin + tol_level = 0.01 + + # A sphere has constant curvature. Both K1 and K2 are equivalent. curvature is 1/r + @testset "Triangulated sphere" begin + r = 10 + k_true = 1.0/r + F,V = geosphere(3,r) + K1,K2,U1,U2,H,G = mesh_curvature_polynomial(F,V) + @test isapprox(mean(K1),k_true,atol=tol_level) + @test isapprox(mean(K2),k_true,atol=tol_level) + @test isapprox(mean(H),k_true,atol=tol_level) + @test isapprox(sqrt(abs(mean(G))),k_true,atol=tol_level) + end + + @testset "Quadrangulated sphere" begin + r = 10 + k_true = 1.0/r + F,V = quadsphere(4,r) + K1,K2,U1,U2,H,G = mesh_curvature_polynomial(F,V) + @test isapprox(mean(K1),k_true,atol=tol_level) + @test isapprox(mean(K2),k_true,atol=tol_level) + @test isapprox(mean(H),k_true,atol=tol_level) + @test isapprox(sqrt(abs(mean(G))),k_true,atol=tol_level) + end + + # For this algorithm the cube appears to have homogeneous curvature too that is 1.5/r + @testset "Cube" begin + r = sqrt(3) + k_true = 1.5*(1.0/r) # Note only for the polynomial method + M = cube(r) + F = faces(M) + V = coordinates(M) + K1,K2,U1,U2,H,G = mesh_curvature_polynomial(F,V) + @test isapprox(mean(K1),k_true,atol=tol_level) + @test isapprox(mean(K2),k_true,atol=tol_level) + @test isapprox(mean(H),k_true,atol=tol_level) + @test isapprox(sqrt(abs(mean(G))),k_true,atol=tol_level) + end + + # For a cylinder k1=1/r while k2=0, hence H = 1/(2*R) + @testset "Cylinder" begin + r = 25 + s = r/10 + nc = round(Int64,(2*pi*r)/s) + d = r*2 + Vc = circlepoints(r, nc; dir=:cw) + num_steps = round(Int64,d/s) + num_steps = num_steps + Int64(iseven(num_steps)) + F, V = extrudecurve(Vc, d; s=1, n=[0.0,0.0,1.0], num_steps=num_steps, close_loop=true, face_type=:quad) + K1,K2,U1,U2,H,G = mesh_curvature_polynomial(F,V) + @test isapprox(mean(K1),1.0/r,atol=tol_level) + @test isapprox(mean(K2),0.0,atol=tol_level) + @test isapprox(mean(H),1.0/(2*r),atol=tol_level) + @test isapprox(sqrt(abs(mean(G))),0.0,atol=tol_level) + end + +end + @testset "separate_vertices" begin eps_level = 0.001