Added tests for mesh_curvature_polynomial #20

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

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) 
 

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