Skip to content

Commit

Permalink
Fixed bugs in connectivity and curvature analysis, added squircle fun…
Browse files Browse the repository at this point in the history
…ction and testing
  • Loading branch information
Kevin-Mattheus-Moerman committed Aug 12, 2024
1 parent 95a8054 commit d709019
Show file tree
Hide file tree
Showing 9 changed files with 382 additions and 86 deletions.
8 changes: 6 additions & 2 deletions examples/demo_circlepoints.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,21 +6,25 @@ using LinearAlgebra
## Defining a set of points on a circle using scalar radius
r = 1.0 # Radius
n = 40 # Number of points
V1 = circlepoints(r,n)
V1 = circlepoints(r,n; dir=:acw)

## Applying a radial function
rFun(t) = r + 0.5.*sin(3*t)
V2 = circlepoints(rFun,n)
V2 = circlepoints(rFun,n; dir=:cw)

## Visualization
fig = Figure(size=(1600,800))

ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A circle of " * string(n) * " points")
lines!(ax1,V1, linewidth = 5, color = :blue)
scatter!(ax1,V1,markersize=25,color = :black)
scatter!(ax1,V1[1],markersize=35,color = :red)
scatter!(ax1,V1[2],markersize=35,color = :yellow)

ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Radial function")
lines!(ax2,V2, linewidth = 5, color = :blue)
scatter!(ax2,V2,markersize=25,color = :black)
scatter!(ax2,V2[1],markersize=35,color = :red)
scatter!(ax2,V2[2],markersize=35,color = :yellow)

fig
43 changes: 28 additions & 15 deletions examples/demo_mesh_curvature_polynomial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,11 @@ using Comodo
using GLMakie
using GeometryBasics
using FileIO
using Rotations
using Statistics

# Example geometry
testCase = 7
# Example geometry0
testCase = 4
if testCase == 1
r = 25.25
F,V = geosphere(5,r)
Expand All @@ -26,7 +28,9 @@ elseif testCase==4
Vc = circlepoints(r, nc; dir=:cw)
num_steps = round(Int,d/s)
num_steps = num_steps + Int(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)
F, V = extrudecurve(Vc; extent=d, n=[0.0,0.0,1.0], num_steps=num_steps, close_loop=true, face_type=:tri)
R = rand(RotXYZ)
V = [R* v for v in V]
elseif testCase==5 # Merged STL for single object
# Loading a mesh
fileName_mesh = joinpath(comododir(),"assets","stl","stanford_bunny_low.stl")
Expand All @@ -36,7 +40,7 @@ elseif testCase==5 # Merged STL for single object
F = tofaces(faces(M))
V = topoints(coordinates(M))
F,V,_ = mergevertices(F,V)
F,V = subtri(F,V,2; method = :loop)
# F,V = subtri(F,V,2; method = :Loop)
elseif testCase==6 # Merged STL for single object
# Loading a mesh
fileName_mesh = joinpath(comododir(),"assets","stl","david.stl")
Expand All @@ -54,36 +58,45 @@ elseif testCase==7
# Obtain mesh faces and vertices
F = tofaces(faces(M))
V = topoints(coordinates(M))
elseif testCase == 8
nSub = 4 # Number of refinement steps of the geodesic sphere
r = 2.5 # Sphere radius
F,V = geosphere(nSub,r) # Creating the faces and vertices of a full sphere
VC = simplexcenter(F,V) # Finding triangle centre coordiantes
F = [F[i] for i in findall(map(v-> v[3]>0,VC))] # Remove some faces using z of central coordinates
F,V = remove_unused_vertices(F,V) # Cleanup/remove unused vertices after faces were removed
end

M = GeometryBasics.Mesh(V,F)

K1,K2,U1,U2,H,G = mesh_curvature_polynomial(F,V)
K1,K2,U1,U2,H,G = mesh_curvature_polynomial(F,V; growsteps=2)

## Visualization
s = pointspacingmean(F,V)
cMap = :Spectral
cMap = Makie.Reverse(:Spectral)
f = 0.1

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))
# 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)
hp1 = poly!(ax1,M, color=K1, strokecolor=:white, strokewidth=0.1, shading = FastShading, transparency=false,colormap=cMap,colorrange = maximum(abs.(filter(!isnan,K1))).*f.*(-1,1))
hpn1 = dirplot(ax1,V,U1; color=:black,linewidth=2,scaleval=s,style=:through)
# scatter!(ax1,V,color=K1,colormap=cMap,colorrange = maximum(abs.(K1)).*0.1.*(-1,1),markersize=10);
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))
hp2 = poly!(ax1,M, color=K2, strokecolor=:white, strokewidth=0.1, shading = FastShading, transparency=false,colormap=cMap,colorrange = maximum(abs.(filter(!isnan,K2))).*f.*(-1,1))
hpn2 = dirplot(ax1,V,U2; color=:black,linewidth=2,scaleval=s,style=:through)
Colorbar(fig[1, 4],hp1)
# scatter!(ax1,V,color=K2,colormap=cMap,colorrange = maximum(abs.(K2)).*0.1.*(-1,1),markersize=10);
Colorbar(fig[1, 4],hp2)

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))
Colorbar(fig[2, 2],hp1)
hp3 = poly!(ax1,M, color=H, strokecolor=:white, strokewidth=0.1, shading = FastShading, transparency=false,colormap=cMap,colorrange = maximum(abs.(filter(!isnan,H))).*f.*(-1,1))
Colorbar(fig[2, 2],hp3)

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))
Colorbar(fig[2, 4],hp1)
hp4 = poly!(ax1,M, color=G, strokecolor=:white, strokewidth=0.1, shading = FastShading, transparency=false,colormap=cMap,colorrange = maximum(abs.(filter(!isnan,G))).*f.*(-1,1))
Colorbar(fig[2, 4],hp4)

fig

Expand Down
13 changes: 4 additions & 9 deletions examples/demo_meshconnectivity.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ using Comodo
using GeometryBasics
using GLMakie

testCase = 1
testCase = 2
if testCase==1
r = 1
M = platonicsolid(4,r) # Get icosahedron
Expand Down Expand Up @@ -37,63 +37,58 @@ con_V2E = C.vertex_edge
con_V2F = C.vertex_face
con_V2V = C.vertex_vertex

indE = 41
indF = 5
indV = 2

fig = Figure(size=(1200,800))

ax=Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "edge-face")
hp=poly!(ax,GeometryBasics.Mesh(V,F), strokewidth=3,color=:white, shading=FastShading, overdraw=false)

indE = 41
hp1=wireframe!(ax,GeometryBasics.Mesh(V,[E_uni[indE]]),linewidth=5, transparency=true, depth_shift=-1.0f-3, color=:blue)
hp2=poly!(ax,GeometryBasics.Mesh(V,F[con_E2F[indE]]), strokewidth=4,color=:red, shading=FastShading, overdraw=false)

ax=Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face-edge")
hp=poly!(ax,GeometryBasics.Mesh(V,F), strokewidth=3,color=:white, shading=FastShading, overdraw=false)

indF = 5
hp1=poly!(ax,GeometryBasics.Mesh(V,[F[indF]]), strokewidth=4,color=:blue, shading=FastShading, overdraw=false)
hp2=wireframe!(ax,GeometryBasics.Mesh(V,E_uni[con_F2E[indF]]),linewidth=5, transparency=true, depth_shift=-1.0f-3, color=:red)

ax=Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face-face (via edges)")
hp=poly!(ax,GeometryBasics.Mesh(V,F), strokewidth=3,color=:white, shading=FastShading, overdraw=false)

indF = 5
hp1=poly!(ax,GeometryBasics.Mesh(V,[F[indF]]), strokewidth=4,color=:blue, shading=FastShading, overdraw=false)
hp2=poly!(ax,GeometryBasics.Mesh(V,F[con_F2F[indF]]), strokewidth=4,color=:red, shading=FastShading, overdraw=false)

ax=Axis3(fig[1, 4], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "vertex-face")
hp=poly!(ax,GeometryBasics.Mesh(V,F), strokewidth=3,color=:white, shading=FastShading, overdraw=false)

indV = 2
hp1=scatter!(ax,V[indV],color =:blue, markersize = 25)
hp2=poly!(ax,GeometryBasics.Mesh(V,F[con_V2F[indV]]), strokewidth=4,color=:red, shading=FastShading, overdraw=false)

ax=Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "vertex-edge")
hp=poly!(ax,GeometryBasics.Mesh(V,F), strokewidth=3,color=:white, shading=FastShading, overdraw=false)

indV = 2
hp1=scatter!(ax,V[indV],color =:blue, markersize = 25)
hp2=wireframe!(ax,GeometryBasics.Mesh(V,E_uni[con_V2E[indV]]),linewidth=5, transparency=true, depth_shift=-1.0f-3, color=:red)


ax=Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "vertex-vertex")
hp=poly!(ax,GeometryBasics.Mesh(V,F), strokewidth=3,color=:white, shading=FastShading, overdraw=false)

indV = 2
hp1=scatter!(ax,V[indV],color =:blue, markersize = 25)
hp2=scatter!(ax,V[con_V2V[indV]],color =:red, markersize = 25)

ax=Axis3(fig[2, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "edge-edge")
hp=poly!(ax,GeometryBasics.Mesh(V,F), strokewidth=3,color=:white, shading=FastShading, overdraw=false)

indE = 41
hp1=wireframe!(ax,GeometryBasics.Mesh(V,[E_uni[indE]]), linewidth=5, transparency=true, depth_shift=-1.0f-3, color=:blue)
hp2=wireframe!(ax,GeometryBasics.Mesh(V,E_uni[con_E2E[indE]]),linewidth=5, transparency=true, depth_shift=-1.0f-3, color=:red)

ax=Axis3(fig[2, 4], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face-face (via vertices)")
hp=poly!(ax,GeometryBasics.Mesh(V,F), strokewidth=3,color=:white, shading=FastShading, overdraw=false)

indF = 5
hp1=poly!(ax,GeometryBasics.Mesh(V,[F[indF]]), strokewidth=4,color=:blue, shading=FastShading, overdraw=false)
hp2=poly!(ax,GeometryBasics.Mesh(V,F[con_F2F_v[indF]]), strokewidth=4,color=:red, shading=FastShading, overdraw=false)

Expand Down
47 changes: 47 additions & 0 deletions examples/demo_squircle.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@
using Comodo
using GeometryBasics
using GLMakie
using Printf

τ = 0.5
n = 4*200
r = 1.0
a_tol = 1e-6



V = squircle(r,n,τ; atol=a_tol)

# Visualization

markersize = 6
linewidth = 10
nSlider = 150
cMap = cgrad(:Spectral, nSlider, categorical = true)#Makie.Reverse(:Spectral)

fig = Figure(size = (1600,1600))
# ax = Axis3(fig[1, 1],aspect = :data,title="Squircle, τ="*string(τ))
ax = Axis(fig[1, 1],aspect = DataAspect(),title="Squircle, τ="*@sprintf("%.3f",τ),titlesize=50)

stepRange1 = range(0.0,1.0,nSlider)

for (i,stepVal) in enumerate(stepRange1)
V = squircle(r,n,stepVal; atol=a_tol)
lines!(ax, V,linewidth=0.25,color=cMap[i],transparency=true,depth_shift=0.01f0)
end

hSlider1 = Slider(fig[2, 1], range = stepRange1, startvalue = 0,linewidth=30)

hp1 = lines!(ax, V,linewidth=linewidth,color=:black,depth_shift=-0.01f0)

on(hSlider1.value) do stepVal
hp1[1] = squircle(r,n,stepVal; atol=a_tol)
ax.title = "Squircle, τ="*@sprintf("%.3f",stepVal)
end

Colorbar(fig[1, 2], limits = (0, 1), colormap = cMap)

fig

# fileName = comododir()*"/assets/img/squircle.mp4"
# slider2anim(fig,hSlider1,fileName; backforth=true, duration=3)
2 changes: 1 addition & 1 deletion examples/demo_sweeploft.jl
Original file line number Diff line number Diff line change
Expand Up @@ -169,7 +169,7 @@ elseif testCase == 4
R = RotMatrix3{Float64}(S22\S21)
V2 = [R*v for v V2] # Rotate the coordinates
V2 = [v .+ Vc[end] for v V2]
elseif testCase ==5
elseif testCase == 5
nc = 75 # Number of points on guide curve
r = 8
a = 1.5*π
Expand Down
68 changes: 68 additions & 0 deletions examples/demo_sweeploft_squircle.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
using Comodo
using GeometryBasics
using GLMakie
using Rotations
using Statistics
using LinearAlgebra

function makeMesh= 1.0)

nc = 150 # Number of points on guide curve
r = 5
Vc = [GeometryBasics.Point{3, Float64}(r*cos(t),r*sin(t),0.0) for t range(0,2*π,nc)]

# Define section curves
np = 40 # Number of section points

# Section 1
V1 = squircle(2.0,np,τ; dir=:acw)
# V1 = evenly_sample(V1, np; close_loop=true)
V1 = [Point{3,Float64}(v[1],v[3],v[2]) for v V1]
V1 = [v .+ Vc[end] for v V1]

# Section 2
V2 = squircle(2.0,np,τ; dir=:acw)
# V2 = evenly_sample(V2, np; close_loop=true)
V2 = [Point{3,Float64}(v[1],v[3],v[2]) for v V2]
V2 = [v .+ Vc[end] for v V2]

F,V = sweeploft(Vc,V1,V2; face_type=:quad, num_twist=3)
F,V = mergevertices(F,V)
# F = invert_faces(F)
# C = ceil.(collect(1:length(V))./(length(V1)))
K1,K2,U1,U2,H,G = mesh_curvature_polynomial(F,V; growsteps=2)

return F,V,G
end

F,V,C = makeMesh(1.0)

#########


# Visualization

cMap = :Spectral

markersize = 6
linewidth = 2

fig = Figure(size = (1600,1600))
ax = Axis3(fig[1, 1],aspect = :data,title="Swept lofting")

stepRange1 = range(0.0,1.0,100)
hSlider1 = Slider(fig[2, 1], range = stepRange1, startvalue = 0,linewidth=30)

# hp1 = poly!(ax, GeometryBasics.Mesh(V,F), strokecolor=:black, strokewidth=0.5,color=C,transparency=false,shading = FastShading,colormap=cMap)
hp1 = mesh!(ax, GeometryBasics.Mesh(V,F), color=C,transparency=false,shading = FastShading,colormap=cMap)

on(hSlider1.value) do τ
F,V,C = makeMesh(τ)
hp1[1] = GeometryBasics.Mesh(V,F)
hp1.color = C
end

fig

# fileName = comododir()*"/assets/img/sweep_loft_squircle_anim_02.mp4"
# slider2anim(fig,hSlider1,fileName; backforth=true, duration=4)
1 change: 1 addition & 0 deletions src/Comodo.jl
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@ export invert_faces, kabsch_rot, sweeploft, loftpoints2surf, revolvecurve
export batman, tridisc, triplate, regiontrimesh, scalesimplex, subcurve, dualclad
export tet2hex, element2faces, subhex, rhombicdodecahedron, tri2quad
export tetgenmesh, surfacevolume, tetvolume, extrudefaces, filletcurve
export squircle, circlerange

# Export functions: Visualization related
export slidercontrol, slider2anim, dirplot, normalplot
Expand Down
Loading

0 comments on commit d709019

Please sign in to comment.