Testing of smoothing methods and quadplate, work towards #20

examples/demo_quadplate.jl

@@ -2,8 +2,8 @@ using Comodo
 using GLMakie
 using GeometryBasics
 
-plateDim = (20,20)
-plateElem = (10,10)
+plateDim = [20,20]
+plateElem = [10,10]
 
 F,V = quadplate(plateDim,plateElem)
 

examples/demo_surface_mesh_smoothing.jl

@@ -20,9 +20,6 @@ V0 = deepcopy(V)
 
 V = map(x-> x.+ (5.0 .* randn(eltype(V))),V)
 
-E_uni = meshedges(F;unique_only=true)
-con_V2V = con_vertex_vertex(E_uni)
-
 # Set smoothing parameters
 tol = 1e-3
 nMax = 100 # Maximum number of iterations
@@ -33,12 +30,20 @@ nMax = 100 # Maximum number of iterations
 α = 0.1
 β = 0.5
 
+testMode = 1
+if testMode == 2
+    z = [v[3] for v in V]
+    constrained_points = findall(z.<0)
+else
+    constrained_points = nothing
+end
+
 ## VISUALISATION
 
 strokeWidth1 = 1
 
-Vs_HC = smoothmesh_hc(F,V,con_V2V; n=nMax,α=α,β=β,tolDist=tol)
-Vs_LAP = smoothmesh_laplacian(F,V,con_V2V; n=nMax,λ=λ)
+Vs_HC = smoothmesh_hc(F,V,nMax,α,β; tolDist=tol, constrained_points = constrained_points)
+Vs_LAP = smoothmesh_laplacian(F,V,nMax,λ; constrained_points = constrained_points)
 Ds = [sqrt(sum((Vs_HC[i]-V0[i]).^2)) for i ∈ eachindex(V)]
 Ds_LAP = [sqrt(sum((Vs_LAP[i]-V0[i]).^2)) for i ∈ eachindex(V)]
 cLim = maximum(Ds_LAP).*(0.0,1.0)
@@ -66,14 +71,14 @@ slidercontrol(hSlider,fig)
 
 on(hSlider.value) do stepIndex
     # Update first plot
-    Vs_LAP = smoothmesh_laplacian(F,V,con_V2V; n=stepIndex,λ=λ)
+    Vs_LAP = smoothmesh_laplacian(F,V,stepIndex,λ; constrained_points = constrained_points)
     Ds_LAP = [sqrt(sum((Vs_LAP[i]-V0[i]).^2)) for i ∈ eachindex(V)]
     ax3.title= "Laplacian smoothed n = " * string(stepIndex) * " times"
     hp3.color = Ds_LAP
     hp3[1] = GeometryBasics.Mesh(Vs_LAP,F)
 
     # Update second plot
-    Vs_HC = smoothmesh_hc(F,V,con_V2V; n=stepIndex,α=α,β=β,tolDist=tol)
+    Vs_HC = smoothmesh_hc(F,V,stepIndex,α,β; con_V2V= con_V2V, tolDist=tol, constrained_points = constrained_points)
     Ds_HC = [sqrt(sum((Vs_HC[i]-V0[i]).^2)) for i ∈ eachindex(V)]
     ax4.title= "HC smoothed n = " * string(stepIndex) * " times"
     hp4.color = Ds_HC

src/functions.jl

@@ -1522,7 +1522,7 @@ function con_vertex_vertex(E,V=nothing,con_V2E=nothing)
     return con_V2V
 end
 
-function meshconnectivity(F,V) #conType = ["ev","ef","ef","fv","fe","ff","ve","vf","vv"]
+function meshconnectivity(F,V)
 
     # EDGE-VERTEX connectivity
     E = meshedges(F)
@@ -1602,14 +1602,14 @@ mean coordinates of that point's Laplacian umbrella. The update features a lerp
 like weighting between the previous iterations coordinates and the mean 
 coordinates. The code features `Vs[q] = (1.0-λ).*Vs[q] .+ λ*mean(V[con_V2V[q]])`
 As can be seen, the weighting is controlled by the input parameter `λ` which is
-in the range (0,1). If `λ=0` then no smoothing occurs. If `λ=1` pure Laplacian mean based smoothing occurs. For 
-intermediate values a linear blending between the two occurs.  
+in the range (0,1). If `λ=0` then no smoothing occurs. If `λ=1` then pure 
+Laplacian mean based smoothing occurs. For intermediate values a linear blending
+between the two occurs.  
 """
-function smoothmesh_laplacian(F,V,con_V2V=nothing; n=1, λ=0.5)
+function smoothmesh_laplacian(F,V,n=1, λ=0.5; con_V2V=nothing, constrained_points=nothing)
     if λ>1.0 || λ<0.0
         error("λ should be in the range 0-1")
     end
-
     if λ>0.0
         # Compute vertex-vertex connectivity i.e. "Laplacian umbrellas" if nothing
         if isnothing(con_V2V)
@@ -1621,6 +1621,9 @@ function smoothmesh_laplacian(F,V,con_V2V=nothing; n=1, λ=0.5)
             for q ∈ eachindex(V)                
                 Vs[q] = (1.0-λ).*Vs[q] .+ λ*mean(V[con_V2V[q]]) # Linear blend between original and pure Laplacian
             end
+            if !isnothing(constrained_points)
+                Vs[constrained_points] = V[constrained_points] # Put back constrained points
+            end
             V = Vs
         end
     end
@@ -1640,7 +1643,7 @@ Laplacian like smoothing but aims to compensate for shrinkage/swelling by also
 # Reference 
 [Vollmer et al. Improved Laplacian Smoothing of Noisy Surface Meshes, 1999. doi: 10.1111/1467-8659.00334](https://doi.org/10.1111/1467-8659.00334)
 """
-function smoothmesh_hc(F,V, con_V2V=nothing; n=1, α=0.1, β=0.5, tolDist=nothing)
+function smoothmesh_hc(F,V, n=1, α=0.1, β=0.5; con_V2V=nothing, tolDist=nothing, constrained_points=nothing)
 
     # Compute vertex-vertex connectivity i.e. "Laplacian umbrellas" if nothing
     if isnothing(con_V2V)
@@ -1658,13 +1661,13 @@ function smoothmesh_hc(F,V, con_V2V=nothing; n=1, α=0.1, β=0.5, tolDist=nothin
             P[i] = mean(Q[con_V2V[i]]) # Laplacian 
             # Compute different vector between P and a point between original 
             # point and Q (which is P before laplacian)
-            B[i] = P[i] .- (α.*V[i] .+ (1-α).*Q[i])
+            B[i] = P[i] .- (α.*V[i] .+ (1.0-α).*Q[i])
         end
         d = 0.0        
         for i ∈ eachindex(V)      
             # Push points back based on blending between pure difference vector
             # B and the Laplacian mean of these      
-            P[i] = P[i] .- (β.*B[i] .+ (1-β).* mean(B[con_V2V[i]]))
+            P[i] = P[i] .- (β.*B[i] .+ (1.0-β).* mean(B[con_V2V[i]]))
         end   
         c+=1 
         if !isnothing(tolDist) # Include tolerance based termination
@@ -1676,6 +1679,9 @@ function smoothmesh_hc(F,V, con_V2V=nothing; n=1, α=0.1, β=0.5, tolDist=nothin
                 break
             end            
         end
+        if !isnothing(constrained_points)
+            P[constrained_points] = V[constrained_points] # Put back constrained points
+        end
     end
     return P
 end

test/runtests.jl

@@ -1048,23 +1048,7 @@ end
 
 end
 
-@testset "mergevertices" begin
-    eps_level = 0.001
-    r = 2 * sqrt(3) / 2
-    M = cube(r)
 
-    F = faces(M)
-    V = coordinates(M)
-    F, V, _ = mergevertices(F, V)
-
-    @test V isa Vector{Point3{Float64}}
-    @test length(V) == 8
-    @test isapprox(V[1], [-1.0, -1.0, -1.0], atol=eps_level)
-
-    @test F isa Vector{QuadFace{Int64}}
-    @test length(F) == 6
-    @test F[1] == [1, 2, 3, 4]
-end
 
 
 @testset "remove_unused_vertices" begin
@@ -1228,6 +1212,75 @@ end
     end
 end
 
+@testset "mergevertices" begin
+    eps_level = 0.001
+    r = 2 * sqrt(3) / 2    
+    M = cube(r)
+    F = faces(M)
+    V = coordinates(M)
+    Fs,Vs = separate_vertices(F,V)
+    Fm, Vm, indReverse = mergevertices(Fs, Vs)
+
+    @test V isa Vector{Point3{Float64}}
+    @test length(Vm) == length(V)
+    @test isapprox(Vm[1], [-1.0, -1.0, -1.0], atol=eps_level)
+    @test [indReverse[f] for f ∈ Fm] == Fs # Reverse mapping 
+    @test Fm isa Vector{QuadFace{Int64}} # Face type unaltered 
+    @test length(Fm) == length(F) # Length is correct
+    @test Fm[1] == [1, 2, 3, 4]
+end
+
+
+@testset "smoothmesh_laplacian" begin
+    eps_level = 0.001
+    M = tetrahedron(1.0)
+    F = faces(M)
+    V = coordinates(M)
+    F,V = subtri(F,V,3)
+
+    λ = 0.5 # Laplacian smoothing parameter
+    n = 10 # Number of iterations 
+    Vs = smoothmesh_laplacian(F,V,n,λ; constrained_points = [5])
+
+    @test Vs[5] == V[5]
+    @test length(V) == length(Vs)
+    @test isapprox(Vs[1:12:end],Point3{Float64}[[-0.5267934833030736, -0.3045704088157244, -0.21536380142235773], [0.18035752833432903, 0.10412811672612346, 0.294521655293559], [0.08303921331720784, 0.4644517405905486, -0.21506879101783555], [0.27578186198836685, 0.061145348885743724, 0.14725778543071683], [0.42141942792006937, -0.14533559613655794, -0.028420418214732533], [0.14563756593170252, 0.013984084814994118, 0.4219980296556624], [0.060704371251411274, 0.5216603792732742, -0.14726169138940381], [0.09542433365403777, 0.33344751143983137, 0.11891253326014106], [0.35882107530557467, -0.013681340938917872, -0.21539751418386643], [0.09542433365403777, 0.22326098199503258, 0.27473981759179783], [0.022334842065796563, 0.5696027951808408, -0.21506313462934012]],atol=eps_level)
+end
+
+
+@testset "smoothmesh_hc" begin
+    eps_level = 0.001
+    M = tetrahedron(1.0)
+    F = faces(M)
+    V = coordinates(M)
+    F,V = subtri(F,V,3)
+
+    n=10
+    α=0.1
+    β=0.5
+    Vs = smoothmesh_hc(F,V,n,α,β; constrained_points = [5])
+
+    @test Vs[5] == V[5]
+    @test length(V) == length(Vs)
+    @test isapprox(Vs[1:12:end],Point3{Float64}[[-0.717172128187643, -0.4136411138967243, -0.29248843661393087], [0.2172725677452308, 0.12542984177391073, 0.354795754721541], [0.14131349634859056, 0.5833163373526694, -0.2928129682724839], [0.3244337573708121, 0.06012744101331828, 0.1773748669098122], [0.5320437804652989, -0.18001157172411225, -0.007876437657606965], [0.2076100230944868, 0.007251571709351666, 0.5218871813767098], [0.09749864497483732, 0.6706350785643729, -0.1774072797262262], [0.10716118962558129, 0.42533928994594383, 0.16951517881969574], [0.4657472537194027, 0.021884259910259368, -0.2924568169614577], [0.10716118962558129, 0.30160020659192405, 0.3445086686945654], [0.04381485137375324, 0.7522177199770612, -0.29279262066755407]],atol=eps_level)
+end
+
+
+@testset "quadplate" begin
+    eps_level = 0.001
+
+    plateDim = [5,5]
+    plateElem = [4,3]
+    F,V = quadplate(plateDim,plateElem)
+
+    @test V isa Vector{Point3{Float64}}
+    @test length(V) == prod(plateElem.+1)
+    @test isapprox(V, Point3{Float64}[[-2.5, -2.5, 0.0], [-1.25, -2.5, 0.0], [0.0, -2.5, 0.0], [1.25, -2.5, 0.0], [2.5, -2.5, 0.0], [-2.5, -0.8333333333333334, 0.0], [-1.25, -0.8333333333333334, 0.0], [0.0, -0.8333333333333334, 0.0], [1.25, -0.8333333333333334, 0.0], [2.5, -0.8333333333333334, 0.0], [-2.5, 0.8333333333333334, 0.0], [-1.25, 0.8333333333333334, 0.0], [0.0, 0.8333333333333334, 0.0], [1.25, 0.8333333333333334, 0.0], [2.5, 0.8333333333333334, 0.0], [-2.5, 2.5, 0.0], [-1.25, 2.5, 0.0], [0.0, 2.5, 0.0], [1.25, 2.5, 0.0], [2.5, 2.5, 0.0]], atol=eps_level)
+    @test F isa Vector{QuadFace{Int64}}
+    @test length(F) == prod(plateElem)
+    @test F[1] == [1, 2, 7, 6]
+end
+
 
 @testset "quadsphere" begin
     r = 1.0