From 03610efacf099d84acf066c2c60d3cdd66028416 Mon Sep 17 00:00:00 2001 From: Kevin-Mattheus-Moerman Date: Sat, 20 Apr 2024 10:55:51 +0100 Subject: [PATCH] Imporved type specifications, added batman curve, and wip for regiontrimesh --- examples/demo_batman.jl | 32 ++ examples/demo_dirplot.jl | 26 + examples/demo_edgeangles.jl | 2 +- examples/demo_facearea.jl | 14 +- examples/demo_facenormal.jl | 1 + examples/demo_kabsch_rot.jl | 8 +- examples/demo_load_view_obj.jl | 12 +- examples/demo_mergevertices.jl | 6 +- examples/demo_mesh_curvature_polynomial.jl | 12 +- examples/demo_meshgroup.jl | 12 +- examples/demo_quad2tri.jl | 4 +- examples/demo_regiontrimesh.jl | 69 +++ examples/demo_remove_unused_vertices.jl | 10 +- examples/demo_subquad_Catmull_Clark.jl | 2 +- examples/demo_subtri_loop.jl | 2 +- examples/demo_surface_mesh_smoothing.jl | 5 +- examples/demo_sweeploft.jl | 74 ++- examples/demo_trisurfslice.jl | 16 +- examples/demo_vertexnormal.jl | 4 +- src/Comodo.jl | 3 +- src/functions.jl | 593 ++++++++++++--------- test/runtests.jl | 189 ++++--- 22 files changed, 692 insertions(+), 404 deletions(-) create mode 100644 examples/demo_batman.jl create mode 100644 examples/demo_dirplot.jl create mode 100644 examples/demo_regiontrimesh.jl diff --git a/examples/demo_batman.jl b/examples/demo_batman.jl new file mode 100644 index 0000000..524a583 --- /dev/null +++ b/examples/demo_batman.jl @@ -0,0 +1,32 @@ +using Comodo +using GLMakie +using DelaunayTriangulation + +#= +This demo shows the use of the `batman` function to create a Batman logo curve. +The curve is useful for testing surface meshing algorithms since it contains +sharp transitions and pointy features. +=# + +n = 75 +V1 = batman(n) + +n = 76 +V2 = batman(n; symmetric = true,dir=:cw) + +# Visualisation +fig = Figure(size=(1800,600)) + +ax1 = Axis3(fig[1, 1],aspect = :data,title="Standard, anti-clockwise") +hp1 = lines!(ax1, V1,linewidth=3,color=:blue) +hp2 = scatter!(ax1, V1,markersize=8,color=:red) +hp2 = scatter!(ax1, V1[1],markersize=15,color=:yellow) +hp2 = scatter!(ax1, V1[2],markersize=15,color=:orange) + +ax2 = Axis3(fig[1, 2],aspect = :data,title="Symmetric, clockwise") +hp1 = lines!(ax2, V2,linewidth=3,color=:blue) +hp2 = scatter!(ax2, V2,markersize=8,color=:red) +hp2 = scatter!(ax2, V2[1],markersize=15,color=:yellow) +hp2 = scatter!(ax2, V2[2],markersize=15,color=:orange) + +fig diff --git a/examples/demo_dirplot.jl b/examples/demo_dirplot.jl new file mode 100644 index 0000000..a31d9c5 --- /dev/null +++ b/examples/demo_dirplot.jl @@ -0,0 +1,26 @@ +using Comodo +using GLMakie +using GeometryBasics +using FileIO +using Statistics + +#= +This demo shows the use of the `dirplot` function to visualize directional data. +=# + +fig = Figure(size=(1600,800)) + +M = icosahedron() +F = faces(M) +V = coordinates(M) +NV = vertexnormal(F,V; weighting=:area) + +styleSet = (:to,:from,:through) + +for i in eachindex(styleSet) + ax1 = Axis3(fig[1, i], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = string(styleSet[i])) + hp1 = poly!(ax1,M, strokewidth=1,shading=FastShading,color=:white, transparency=true, overdraw=false) + hpa = dirplot(ax1,V,NV.*mean(edgelengths(F,V)); color=:blue,linewidth=3,scaleval=1.0,style=styleSet[i]) +end + +fig diff --git a/examples/demo_edgeangles.jl b/examples/demo_edgeangles.jl index 84bffc1..7acb6d1 100644 --- a/examples/demo_edgeangles.jl +++ b/examples/demo_edgeangles.jl @@ -14,7 +14,7 @@ end a = pi/4 # "45 degree shear" f[1,2] = tan(a) -V2 = togeometrybasics_points([f*v for v ∈ V]) # Shear the cube +V2 = [eltype(V)(f*v) for v ∈ V] # Shear the cube A = edgeangles(F,V) A2 = edgeangles(F,V2) diff --git a/examples/demo_facearea.jl b/examples/demo_facearea.jl index 596f558..af426fd 100644 --- a/examples/demo_facearea.jl +++ b/examples/demo_facearea.jl @@ -31,18 +31,18 @@ elseif testCase==4 # Merged STL for single object M = load(fileName_mesh) # Obtain mesh faces and vertices - F = togeometrybasics_faces(faces(M)) - V = togeometrybasics_points(coordinates(M)) + 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,3; method = :Loop) elseif testCase==5 # Merged STL for single object # Loading a mesh fileName_mesh = joinpath(comododir(),"assets","stl","david.stl") M = load(fileName_mesh) # Obtain mesh faces and vertices - F = togeometrybasics_faces(faces(M)) - V = togeometrybasics_points(coordinates(M)) + F = tofaces(faces(M)) + V = topoints(coordinates(M)) F,V,_ = mergevertices(F,V) elseif testCase==6 # Loading a mesh @@ -50,8 +50,8 @@ elseif testCase==6 M = load(fileName_mesh) # Obtain mesh faces and vertices - F = togeometrybasics_faces(faces(M)) - V = togeometrybasics_points(coordinates(M)) + F = tofaces(faces(M)) + V = topoints(coordinates(M)) end diff --git a/examples/demo_facenormal.jl b/examples/demo_facenormal.jl index 8344a61..7e79411 100644 --- a/examples/demo_facenormal.jl +++ b/examples/demo_facenormal.jl @@ -1,5 +1,6 @@ using Comodo using GLMakie +using GeometryBasics #= This demo shows the use of the `meshnormal` function to obtain mesh face normal diff --git a/examples/demo_kabsch_rot.jl b/examples/demo_kabsch_rot.jl index 1599181..1312f55 100644 --- a/examples/demo_kabsch_rot.jl +++ b/examples/demo_kabsch_rot.jl @@ -11,12 +11,18 @@ are visualized. =# # Loading a mesh + +# fileName_mesh = joinpath(comododir(),"assets","obj","spot_control_mesh_texture.obj") +# M = load(fileName_mesh) + fileName_mesh = joinpath(comododir(),"assets","stl","stanford_bunny_low.stl") M = load(fileName_mesh) # Obtain mesh faces and vertices F = faces(M) -V1 = togeometrybasics_points(coordinates(M)) +V1 = topoints(coordinates(M)) +V1 = [Point{3, Float64}(v) for v in V1] + F,V1 = mergevertices(F,V1) # Define a rotation tensor using Euler angles diff --git a/examples/demo_load_view_obj.jl b/examples/demo_load_view_obj.jl index 9a18ab3..9fcb2e7 100644 --- a/examples/demo_load_view_obj.jl +++ b/examples/demo_load_view_obj.jl @@ -6,22 +6,16 @@ using FileIO fileName_mesh = joinpath(comododir(),"assets","obj","spot_quadrangulated.obj") # M = load(fileName_mesh; facetype = GeometryBasics.QuadFace{GeometryBasics.GLIndex}) M = load(fileName_mesh) -F = togeometrybasics_faces(faces(M)) -V = togeometrybasics_points(coordinates(M)) +F = tofaces(faces(M)) +V = topoints(coordinates(M)) # M = GeometryBasics.Mesh(V,F) -# F,V = mergevertices(F,V) - -# M=cube(1) - - ## Visualization fig = Figure(size=(800,800)) ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Spot the cow") -hp1 = poly!(ax1,M, color=load(joinpath(comododir(),"assets","obj","spot_texture.png")), shading = FastShading, transparency=false,strokecolor=:black,strokewidth=2) +hp1 = poly!(ax1,M, color=load(joinpath(comododir(),"assets","obj","spot_texture.png")), shading = FastShading, transparency=false,strokecolor=:black,strokewidth=1) # hp1 = mesh!(ax1,M, color=:red, shading = FastShading, transparency=false) # hp1 = poly!(ax1,M, color=:red, shading = FastShading, transparency=false,strokecolor=:black,strokewidth=2) # normalplot(ax1,M) - fig \ No newline at end of file diff --git a/examples/demo_mergevertices.jl b/examples/demo_mergevertices.jl index 80cec71..10ed8f3 100644 --- a/examples/demo_mergevertices.jl +++ b/examples/demo_mergevertices.jl @@ -8,10 +8,8 @@ using FileIO # Loading a mesh fileName_mesh = joinpath(comododir(),"assets","stl","stanford_bunny_low.stl") M1 = load(fileName_mesh) -F1 = faces(M1) -V1 = coordinates(M1) -F1 = togeometrybasics_faces(F1) -V1 = togeometrybasics_points(V1) +F1 = tofaces(faces(M1)) +V1 = topoints(coordinates(M1)) F2 = deepcopy(F1) V2 = deepcopy(V1) diff --git a/examples/demo_mesh_curvature_polynomial.jl b/examples/demo_mesh_curvature_polynomial.jl index ae0877c..11ff0ca 100644 --- a/examples/demo_mesh_curvature_polynomial.jl +++ b/examples/demo_mesh_curvature_polynomial.jl @@ -33,8 +33,8 @@ elseif testCase==5 # Merged STL for single object M = load(fileName_mesh) # Obtain mesh faces and vertices - F = togeometrybasics_faces(faces(M)) - V = togeometrybasics_points(coordinates(M)) + F = tofaces(faces(M)) + V = topoints(coordinates(M)) F,V,_ = mergevertices(F,V) F,V = subtri(F,V,2; method = :loop) elseif testCase==6 # Merged STL for single object @@ -43,8 +43,8 @@ elseif testCase==6 # Merged STL for single object M = load(fileName_mesh) # Obtain mesh faces and vertices - F = togeometrybasics_faces(faces(M)) - V = togeometrybasics_points(coordinates(M)) + F = tofaces(faces(M)) + V = topoints(coordinates(M)) F,V,_ = mergevertices(F,V) elseif testCase==7 # Loading a mesh @@ -52,8 +52,8 @@ elseif testCase==7 M = load(fileName_mesh) # Obtain mesh faces and vertices - F = togeometrybasics_faces(faces(M)) - V = togeometrybasics_points(coordinates(M)) + F = tofaces(faces(M)) + V = topoints(coordinates(M)) end M = GeometryBasics.Mesh(V,F) diff --git a/examples/demo_meshgroup.jl b/examples/demo_meshgroup.jl index 14edb02..fcc9956 100644 --- a/examples/demo_meshgroup.jl +++ b/examples/demo_meshgroup.jl @@ -66,8 +66,8 @@ elseif testCase==5 # Unmerged STL, each triangle is seperate group M = load(fileName_mesh) # Obtain mesh faces and vertices - F = togeometrybasics_faces(faces(M)) - V = togeometrybasics_points(coordinates(M)) + F = tofaces(faces(M)) + V = topoints(coordinates(M)) elseif testCase==6 # Merged STL for single object # Loading a mesh using FileIO @@ -75,8 +75,8 @@ elseif testCase==6 # Merged STL for single object M = load(fileName_mesh) # Obtain mesh faces and vertices - F = togeometrybasics_faces(faces(M)) - V = togeometrybasics_points(coordinates(M)) + F = tofaces(faces(M)) + V = topoints(coordinates(M)) F,V,_ = mergevertices(F,V) elseif testCase==7 # Merged STL for single object # Loading a mesh @@ -85,8 +85,8 @@ elseif testCase==7 # Merged STL for single object M = load(fileName_mesh) # Obtain mesh faces and vertices - F = togeometrybasics_faces(faces(M)) - V = togeometrybasics_points(coordinates(M)) + F = tofaces(faces(M)) + V = topoints(coordinates(M)) F,V,_ = mergevertices(F,V) F2 = map(f-> f.+length(V),deepcopy(F)) diff --git a/examples/demo_quad2tri.jl b/examples/demo_quad2tri.jl index fc9dbc1..a3666d5 100644 --- a/examples/demo_quad2tri.jl +++ b/examples/demo_quad2tri.jl @@ -34,7 +34,7 @@ elseif testCase ==3 end a = pi/6 # "45 degree shear" f[1,2] = tan(a) - V = togeometrybasics_points([f*v for v ∈ V]) # Shear the cube + V = [eltype(V)(f*v) for v ∈ V] # Shear the cube # Build deformation gradient tensor to induce shear with known angles f = zeros(3,3) @@ -44,7 +44,7 @@ elseif testCase ==3 a = pi/6 # "45 degree shear" f[3,1] = tan(a) - V = togeometrybasics_points([f*v for v ∈ V]) # Shear the cube + V = [eltype(V)(f*v) for v ∈ V] # Shear the cube end diff --git a/examples/demo_regiontrimesh.jl b/examples/demo_regiontrimesh.jl new file mode 100644 index 0000000..cc7eff7 --- /dev/null +++ b/examples/demo_regiontrimesh.jl @@ -0,0 +1,69 @@ +using Comodo +using GLMakie +using DelaunayTriangulation + +#= +This demo shows the use of `hexbox` to generate a hexahedral mesh for a 3D box +domain. +=# + +n = 100 +V1 = batman(n) + +V2 = deepcopy(V1) +Vc = circlepoints(0.25,25; dir=:cw) +V2 = append!(V2,Vc) +V2 = [Point{2,Float64}(v[1],v[2]) for v in V2] + +TR_1 = triangulate(V1) +F1 = [TriangleFace{Int64}(tr) for tr in TR_1.triangles] + +constrained_segments = append!(collect(1:length(V1)),1) +TR_2 = triangulate(V1; boundary_nodes=constrained_segments) +F2 = [TriangleFace{Int64}(tr) for tr in TR_2.triangles] + +constrained_segments1 = append!(collect(1:length(V1)),1) +constrained_segments2 = append!(collect(1:length(Vc)),1).+length(V1) +constrained_segments = [[constrained_segments1], [constrained_segments2]] +TR_3 = triangulate(V2; boundary_nodes=constrained_segments) +F3 = [TriangleFace{Int64}(tr) for tr in TR_3.triangles] +V3 = [Point{3,Float64}(v[1],v[2],0.0) for v in V2] + +TR_4 = deepcopy(TR_3) +stats = refine!(TR_4; max_area=1e-3, min_angle=27.3) +F4 = [TriangleFace{Int64}(tr) for tr in TR_4.triangles] +V4 = TR_4.points +V4 = [Point{3,Float64}(v[1],v[2],0.0) for v in V4] + +E = boundaryedges(F4) +indB = unique(reduce(vcat,E)) +n = 25 +λ = 0.5 +V4 = smoothmesh_laplacian(F4,V4, n, λ; constrained_points=indB) + +# Visualisation +fig = Figure(size=(1800,600)) + +ax1 = Axis3(fig[1, 1],aspect = :data,title="Unconstrained") +hp1 = lines!(ax1, V1,linewidth=4,color=:blue) +hp2 = scatter!(ax1, V1,markersize=8,color=:red) +hp3 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=0.5,color=:white, strokecolor=:black, shading = FastShading, transparency=false) + +ax2 = Axis3(fig[1, 2],aspect = :data,title="Constrained") +hp1 = lines!(ax2, V1,linewidth=4,color=:blue) +hp2 = scatter!(ax2, V1,markersize=8,color=:red) +hp3 = poly!(ax2,GeometryBasics.Mesh(V1,F2), strokewidth=0.5,color=:white, strokecolor=:black, shading = FastShading, transparency=false) + +ax3 = Axis3(fig[2, 1],aspect = :data,title="Constrained with hole") +hp11 = lines!(ax3, V3[constrained_segments1],linewidth=4,color=:blue) +hp12 = lines!(ax3, V3[constrained_segments2],linewidth=4,color=:blue) +hp2 = scatter!(ax3, V3,markersize=8,color=:red) +hp3 = poly!(ax3,GeometryBasics.Mesh(V3,F3), strokewidth=0.5,color=:white, strokecolor=:black, shading = FastShading, transparency=false) + +ax4 = Axis3(fig[2, 2],aspect = :data,title="Constrained with hole") +hp11 = lines!(ax4, V4[constrained_segments1],linewidth=4,color=:blue) +hp12 = lines!(ax4, V4[constrained_segments2],linewidth=4,color=:blue) +hp2 = scatter!(ax4, V4,markersize=8,color=:red) +hp3 = poly!(ax4,GeometryBasics.Mesh(V4,F4), strokewidth=0.5,color=:white, strokecolor=:black, shading = FastShading, transparency=false) + +fig diff --git a/examples/demo_remove_unused_vertices.jl b/examples/demo_remove_unused_vertices.jl index 17964b8..f01351f 100644 --- a/examples/demo_remove_unused_vertices.jl +++ b/examples/demo_remove_unused_vertices.jl @@ -14,13 +14,15 @@ F = [F[i] for i in findall(map(v-> v[3]>0,VC))] # Remove some faces Fc,Vc = remove_unused_vertices(F,V) - +# Visualization +markersize = 12 M = GeometryBasics.Mesh(Vc,Fc) fig = Figure(size=(1200,1200)) -ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A sliced mesh") - +ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A cut mesh") hp3 = poly!(ax1,M, strokewidth=1,color=:white, strokecolor=:blue, shading = FastShading, transparency=false) -hp3 = normalplot(ax1,M) +# hp3 = normalplot(ax1,M) +scatter!(Vc,markersize=markersize,color=:red) + fig diff --git a/examples/demo_subquad_Catmull_Clark.jl b/examples/demo_subquad_Catmull_Clark.jl index 2029cd6..0c21a23 100644 --- a/examples/demo_subquad_Catmull_Clark.jl +++ b/examples/demo_subquad_Catmull_Clark.jl @@ -20,7 +20,7 @@ rs3 = mean(Rn3) d3 = Rn3.-rs3 ## Visualization -fig = Figure(size=(1600,800)) +fig = Figure(size=(800,800)) ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Refined n=1") wireframe!(ax1,M, linewidth=8,color=:red, overdraw=false) diff --git a/examples/demo_subtri_loop.jl b/examples/demo_subtri_loop.jl index d02893d..75f6778 100644 --- a/examples/demo_subtri_loop.jl +++ b/examples/demo_subtri_loop.jl @@ -27,7 +27,7 @@ ax = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", slidercontrol(hSlider,ax) hp1=wireframe!(ax,M,linewidth=3,color=:red, overdraw=false) -hp2=poly!(ax,Mn,strokewidth=2,color=:white, shading = FastShading) +hp2=poly!(ax,Mn,strokewidth=1,color=:white, shading = FastShading) # hp2=poly!(scene,M,strokewidth=1,color=:white, shading = FastShading) diff --git a/examples/demo_surface_mesh_smoothing.jl b/examples/demo_surface_mesh_smoothing.jl index bd8dc4a..5f1edeb 100644 --- a/examples/demo_surface_mesh_smoothing.jl +++ b/examples/demo_surface_mesh_smoothing.jl @@ -8,9 +8,8 @@ if testCase == 1 # Triangle mesh bunny fileName_mesh = joinpath(comododir(),"assets","stl","stanford_bunny_low.stl") M = load(fileName_mesh) F = faces(M) - V = coordinates(M) - F = togeometrybasics_faces(F) - V = togeometrybasics_points(V) + V = coordinates(M) + V = topoints(V) F,V,ind1,ind2 = mergevertices(F,V; roundVertices=true) elseif testCase == 2 # Quad mesh sphere F,V = quadsphere(4,100) diff --git a/examples/demo_sweeploft.jl b/examples/demo_sweeploft.jl index d10b630..a628d2b 100644 --- a/examples/demo_sweeploft.jl +++ b/examples/demo_sweeploft.jl @@ -9,7 +9,7 @@ testCase = 1 if testCase == 1 # Define guide curve - nc = 51 # Number of points on guide curve + nc = 101 # Number of points on guide curve P = Vector{GeometryBasics.Point{3, Float64}}(undef,4) P[1 ] = GeometryBasics.Point{3, Float64}( 0.0, 0.0, 0.0) P[2 ] = GeometryBasics.Point{3, Float64}( 1.0, 0.0, 0.0) @@ -20,7 +20,7 @@ if testCase == 1 Vc,Sc = evenly_sample(Vc, nc) # Define section curves - np = 25 # Number of section points + np = 50 # Number of section points f(t) = 2.0 + 0.5.*sin(3*t) V1 = circlepoints(f,np; dir=:acw) V1,_ = evenly_sample(V1, np) @@ -110,6 +110,66 @@ elseif testCase == 3 R = RotMatrix3{Float64}(S22\S21) V2 = [R*v for v ∈ V2] # Rotate the coordinates V2 = [v .+ Vc[end] for v ∈ V2] + +elseif testCase == 4 + nc = 1001 # Number of points on guide curve + t = range(0,4.0*pi,nc) + ff = 8.0 + rb = 6.0 + rc = rb/3 + r = [rb + rc*sin(ff*tt) for tt in t] + Vc = [GeometryBasics.Point{3, Float64}(r[i]*cos(t[i]),r[i]*sin(t[i]),t[i]+rc*cos(ff*t[i])) for i in eachindex(t)] + Vc,_ = evenly_sample(Vc, nc) + + # Define section curves + np = 150 # Number of section points + + # Section 1 + f(t) = rc/3 + rc/5 .* sin(3*t) + V1 = circlepoints(f,np; dir=:cw) + V1,_ = evenly_sample(V1, np) + Q = RotXYZ(0.5*π,0.0,0.0) # Define a rotation tensor using Euler angles + V1 = [Q*v for v ∈ V1] # Rotate the coordinates + + # Ensure section is orthogonal to guide curve + n3_1 = Q*Vec3{Float64}(0.0,0.0,-1.0) + n2_1 = Q*Vec3{Float64}(1.0,0.0,0.0) + n1_1 = normalizevector(cross(n3_1,n2_1)) + S11 = mapreduce(permutedims,vcat,[n1_1,n2_1,n3_1]) + + n3_1 = normalizevector(normalizevector(Vc[2]-Vc[1])) + n2_1 = normalizevector(cross(normalizevector(V1[1]),n3_1)) + n1_1 = normalizevector(cross(n3_1,n2_1)) + S12 = mapreduce(permutedims,vcat,[n1_1,n2_1,n3_1]) + + R = RotMatrix3{Float64}(S12\S11) + V1 = [R*v for v ∈ V1] + V1= [v .+ Vc[1] for v ∈ V1] + + # Section 2 + f(t) = rc/4 + rc/6*sin(5*t) + V2 = circlepoints(f,np; dir=:cw) + V2,_ = evenly_sample(V2, np) + Q1 = RotXYZ(0.5*π,0.0,0.0) # Define a rotation tensor using Euler angles + Q2 = RotXYZ(0.0,-0.25*π,0.0) # Define a rotation tensor using Euler angles + Q = Q2*Q1 + V2 = [Q*v for v ∈ V2] # Rotate the coordinates + + # Ensure section is orthogonal to guide curve + n3_2 = Q*Vec3{Float64}(0.0,0.0,-1.0) + n2_2 = Q*Vec3{Float64}(1.0,0.0,0.0) + n1_2 = normalizevector(cross(n3_2,n2_2)) + S21 = RotMatrix3{Float64}(mapreduce(permutedims,vcat,[n1_2,n2_2,n3_2])) + + n3_2 = normalizevector(normalizevector(Vc[2]-Vc[1])) + n2_2 = normalizevector(cross(normalizevector(V1[1]),n3_1)) + n1_2 = normalizevector(cross(n3_1,n2_1)) + S22 = mapreduce(permutedims,vcat,[n1_1,n2_1,n3_1]) + + R = RotMatrix3{Float64}(S22\S21) + V2 = [R*v for v ∈ V2] # Rotate the coordinates + V2 = [v .+ Vc[end] for v ∈ V2] + end @@ -119,13 +179,13 @@ end # face_type=:quad F,V = sweeploft(Vc,V1,V2; face_type=:quad, num_twist=0) F = invert_faces(F) - +C = ceil.(collect(1:length(V))./length(Vc)) # Visualization fig = Figure(size = (800,800)) ax = Axis3(fig[1, 1],aspect = :data,title="Swept lofting") -stepRange1 = -4:1:4 +stepRange1 = -10:1:10 hSlider1 = Slider(fig[2, 1], range = stepRange1, startvalue = 0,linewidth=30) # scatter!(ax, Vc,markersize=8,color=:black) @@ -137,13 +197,13 @@ hp2 = lines!(ax, V1,linewidth=4,color=:blue) scatter!(ax, V2,markersize=8,color=:red) hp3 = lines!(ax, V2,linewidth=4,color=:red) -hp1 = poly!(ax, GeometryBasics.Mesh(V,F), strokecolor=:black, strokewidth=1,color=:white,transparency=false,shading = FastShading) -# hp1 = mesh!(ax, GeometryBasics.Mesh(V,F), color=:white,transparency=false,shading = FastShading) +# hp1 = poly!(ax, GeometryBasics.Mesh(V,F), strokecolor=:black, strokewidth=0.5,color=C,transparency=false,shading = FastShading,colormap=:Spectral) +hp1 = mesh!(ax, GeometryBasics.Mesh(V,F), color=C,transparency=false,shading = FastShading,colormap=Makie.Reverse(:Spectral)) on(hSlider1.value) do stepIndex1 F,V = sweeploft(Vc,V1,V2; face_type=:quad, num_twist=stepIndex1) F = invert_faces(F) - hp1[1] = GeometryBasics.Mesh(V,F) + hp1[1] = GeometryBasics.Mesh(V,F) end fig diff --git a/examples/demo_trisurfslice.jl b/examples/demo_trisurfslice.jl index b731d4e..88853f1 100644 --- a/examples/demo_trisurfslice.jl +++ b/examples/demo_trisurfslice.jl @@ -15,10 +15,8 @@ elseif testCase == 2 M = load(fileName_mesh) # Obtain mesh faces and vertices - F = faces(M) - V = coordinates(M) - F = togeometrybasics_faces(F) - V = togeometrybasics_points(V) + F = tofaces(faces(M)) + V = topoints(coordinates(M)) F,V,ind1,ind2 = mergevertices(F,V) # F,V=subtri(F,V,1) elseif testCase == 3 @@ -27,10 +25,8 @@ elseif testCase == 3 M = load(fileName_mesh) # Obtain mesh faces and vertices - F = faces(M) - V = coordinates(M) - F = togeometrybasics_faces(F) - V = togeometrybasics_points(V) + F = tofaces(faces(M)) + V = topoints(coordinates(M)) F,V,ind1,ind2 = mergevertices(F,V) end @@ -61,7 +57,7 @@ MG = GeometryBasics.Mesh(VGn,FG1) fig = Figure(size=(800,800)) ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A sliced mesh") -stepRange = range(-s,s,50) +stepRange = range(-s,s,100) hSlider = Slider(fig[2, 1], range = stepRange, startvalue = 0,linewidth=30) # hp1 = mesh!(ax1,GeometryBasics.Mesh(V,F),color=:white, shading = FastShading, transparency=true) @@ -84,7 +80,7 @@ on(hSlider.value) do stepIndex end VGn = [GeometryBasics.Point{3, Float64}(R'*v)+pp for v ∈ VG1] # Rotate plane - MG = GeometryBasics.Mesh(togeometrybasics_points(VGn),FG1) + MG = GeometryBasics.Mesh(VGn,FG1) hp2[1] = MG hp3[1] = Mn diff --git a/examples/demo_vertexnormal.jl b/examples/demo_vertexnormal.jl index b1a9475..9cfd556 100644 --- a/examples/demo_vertexnormal.jl +++ b/examples/demo_vertexnormal.jl @@ -34,8 +34,8 @@ for q=1:1:4 M = load(fileName_mesh) # Obtain mesh faces and vertices - F = togeometrybasics_faces(faces(M)) - V = togeometrybasics_points(coordinates(M)) + F = faces(M) + V = topoints(coordinates(M)) F,V,_ = mergevertices(F,V) M = GeometryBasics.Mesh(V,F) titleString="bunny" diff --git a/src/Comodo.jl b/src/Comodo.jl index 8a6694d..199d83a 100644 --- a/src/Comodo.jl +++ b/src/Comodo.jl @@ -20,7 +20,7 @@ export gridpoints, interp_biharmonic_spline, interp_biharmonic, nbezier, gunique export ind2sub, sub2ind, meshedges, edgelengths, unique_simplices, unique_dict export subtri, subquad, geosphere, quad2tri export icosahedron, tetrahedron, cube, dodecahedron, octahedron, platonicsolid -export togeometrybasics_faces, togeometrybasics_points, togeometrybasics_mesh +export tofaces, topoints, togeometrybasics_mesh export quadplate, quadsphere, smoothmesh_laplacian, smoothmesh_hc export mergevertices, separate_vertices, remove_unused_vertices export edgecrossproduct, facearea, facenormal, vertexnormal, normalizevector, dirplot, normalplot @@ -32,6 +32,7 @@ export wrapindex, edgeangles, count_edge_face, boundaryedges, edges2curve export pointspacingmean, extrudecurve, meshgroup, ray_triangle_intersect export distmarch, mesh_curvature_polynomial, curve_length, evenly_sample export invert_faces, kabsch_rot, sweeploft, loftpoints2surf, revolvecurve +export batman # Export functions: Visualization related export slidercontrol diff --git a/src/functions.jl b/src/functions.jl index 62b8ed2..87aeecb 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -119,10 +119,10 @@ end The `gridpoints` function returns a vector of 3D points which span a grid in 3D space. Points are defined as per the input ranges or range vectors. The output -point vector contains elements of the type `GeometryBasics.Point3`. +point vector contains elements of the type `Point`. """ function gridpoints(x::Union{Vector{T}, AbstractRange{T}}, y=x, z=x) where T<:Real - reshape([GeometryBasics.Point{3, Float64}(x, y, z) for z in z, y in y, x in x], + reshape([Point{3, Float64}(x, y, z) for z in z, y in y, x in x], length(x)*length(y)*length(z)) end @@ -131,7 +131,7 @@ end # Description -This function uses biharmonic spline interpolation, which features radial basis +This function uses biharmonic spline interpolation [1], which features radial basis functions. The input is assumed to represent ordered data, i.e. consequtive unique points on a curve. The curve x-, and y-coordinates are provided through the input parameters `x` and `y` respectively. The third input `xi` defines the @@ -139,7 +139,7 @@ sites at which to interpolate. Each of in the input parameters can be either a vector or a range. # References -[David T. Sandwell, Biharmonic spline interpolation of GEOS-3 and SEASAT altimeter data, Geophysical Research Letters, 2, 139-142, 1987. doi: 10.1029/GL014i002p00139](https://doi.org/10.1029/GL014i002p00139) +1. [David T. Sandwell, Biharmonic spline interpolation of GEOS-3 and SEASAT altimeter data, Geophysical Research Letters, 2, 139-142, 1987. doi: 10.1029/GL014i002p00139](https://doi.org/10.1029/GL014i002p00139) """ function interp_biharmonic_spline(x::Union{Vector{T}, AbstractRange{T}},y::Union{Vector{T}, AbstractRange{T}},xi::Union{Vector{T}, AbstractRange{T}}; extrapolate_method=:linear,pad_data=:linear) where T<:Real @@ -236,13 +236,12 @@ end # Description -This function uses biharmonic interpolation. The input `x` should define a +This function uses biharmonic interpolation [1]. The input `x` should define a vector consisting of m points which are n-dimensional, and the input `y` should be a vector consisting of m scalar data values. # References - -[David T. Sandwell, Biharmonic spline interpolation of GEOS-3 and SEASAT altimeter data, Geophysical Research Letters, 2, 139-142, 1987. doi: 10.1029/GL014i002p00139](https://doi.org/10.1029/GL014i002p00139) +1. [David T. Sandwell, Biharmonic spline interpolation of GEOS-3 and SEASAT altimeter data, Geophysical Research Letters, 2, 139-142, 1987. doi: 10.1029/GL014i002p00139](https://doi.org/10.1029/GL014i002p00139) """ function interp_biharmonic(x,y,xi) # Distances from all points in X to all points in X @@ -269,9 +268,9 @@ end This function returns `n` points for an m-th order Bézier spline, based on the m control points contained in the input vector `P`. This function supports point vectors with elements of the type `AbstractPoint{3}` (e.g. -`GeometryBasics.Point{3, Float64}`) or `Vector{Float64}`. +`Point{3, Float64}`) or `Vector{Float64}`. """ -function nbezier(P::Vector{T},n::Integer) where T<:Union{AbstractPoint{3},Vector{Float64}} +function nbezier(P::Vector{Point{ND,TV}},n::Integer) where ND where TV<:Real if n<2 error("The vale of n is too low. Request at least two data points") end @@ -281,13 +280,8 @@ function nbezier(P::Vector{T},n::Integer) where T<:Union{AbstractPoint{3},Vector nnr = N-1:-1:0 f = factorial.(nn) s = factorial(N-1)./( f.*reverse(f) ); # Sigma - - if T<:AbstractPoint{3} - V = [T(0.0,0.0,0.0) for _ in 1:n] - else - V = fill(zeros(Float64,3),n) - end - + + V = zeros(Point{ND,TV},n) @inbounds for i in 1:n b = s.* ((1.0-t[i]).^nnr) .* (t[i].^nn) @inbounds for j = 1:N @@ -345,7 +339,7 @@ end dist(V1,V2) # Description Function compute an nxm distance matrix for the n inputs points in `V1`, and the -m input points in `V2`. The input points may be multidimensional, in face they +m input points in `V2`. The input points may be multidimensional, in fact they can be any type supported by the `euclidean` function of `Distances.jl`. See also: https://github.com/JuliaStats/Distances.jl """ @@ -812,31 +806,24 @@ end """ meshedges(S; unique_only=false) This function returns the edges `E` for the input "simplices" (e.g. faces) -defined by `S`. The input `S` can either represent a vector of faces or a +defined by `F`. The input `F` can either represent a vector of faces or a GeometryBasics mesh. """ -function meshedges(S; unique_only=false) - # Get number of nodes per simplex, use first for now (better to get this from some property) - s1 = S[1] # First simplex - m = length(s1) #Number of nodes per simplex - - E = GeometryBasics.LineFace{Int}[] - for j1 in 1:m # Loop over each node/point for the current simplex - if j11.0 || λ<0.0 throw(ArgumentError("λ should be in the range 0-1")) end @@ -1778,7 +1788,7 @@ function smoothmesh_laplacian(F,V,n=1, λ=0.5; con_V2V=nothing, constrained_poin V = Vs end else #n<0 - throw(ArgumentError("n should greater or equal to 0")) + throw(ArgumentError("n should be greater or equal to 0")) end end return V @@ -1790,14 +1800,14 @@ end # Description -This function implements HC (Humphrey's Classes) smoothing. This method uses +This function implements HC (Humphrey's Classes) smoothing [1]. This method uses Laplacian like smoothing but aims to compensate for shrinkage/swelling by also "pushing back" towards the original coordinates. # 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) +1. [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, n=1, α=0.1, β=0.5; con_V2V=nothing, tolDist=nothing, constrained_points=nothing) +function smoothmesh_hc(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}, n=1, α=0.1, β=0.5; con_V2V=nothing, tolDist=nothing, constrained_points=nothing) where N where TF<:Integer where ND where TV<:Real if α>1.0 || α<0.0 throw(ArgumentError("α should be in the range 0-1")) @@ -1859,10 +1869,10 @@ end function quadplate(plateDim,plateElem) num_x = plateElem[1]+1 num_y = plateElem[2]+1 - V = Vector{GeometryBasics.Point{3, Float64}}() + V = Vector{Point{3, Float64}}() for y = range(-plateDim[2]/2,plateDim[2]/2,num_y) for x = range(-plateDim[1]/2,plateDim[1]/2,num_x) - push!(V,GeometryBasics.Point{3, Float64}(x,y,0.0)) + push!(V,Point{3, Float64}(x,y,0.0)) end end @@ -1877,7 +1887,7 @@ function quadplate(plateDim,plateElem) return F, V end -function quadsphere(n,r) +function quadsphere(n::Int64,r::T) where T <: Real M = platonicsolid(2,r) F = faces(M) V = coordinates(M) @@ -1924,33 +1934,33 @@ function vertex2simplexdata(F,DV) return DF end -function simplexcenter(F,V) +function simplexcenter(F,V::Vector{Point{ND,TV}}) where ND where TV<:Real return vertex2simplexdata(F,V) end -function normalizevector(A) - if eltype(A) <: Number - return A./norm(A) - else +function normalizevector(A::Union{Vector{Point{ND,TV}},Vector{Vec{ND,TV}}}) where ND where TV<:Real return A./norm.(A) - end end -function circlepoints(r,n; dir=:acw) +function normalizevector(A::Union{Point{ND,TV},Vec{ND,TV}}) where ND where TV<:Real + return A./norm(A) +end + +function circlepoints(r::T,n::Int64; dir=:acw) where T <: Real if dir==:acw - return [GeometryBasics.Point{3, Float64}(r*cos(t),r*sin(t),0) for t in range(0.0,2.0*π-(2.0*π)/n,n)] + return [Point{3, Float64}(r*cos(t),r*sin(t),0) for t in range(0.0,2.0*π-(2.0*π)/n,n)] elseif dir==:cw - return [GeometryBasics.Point{3, Float64}(r*cos(t),r*sin(t),0) for t in range(0.0,(2.0*π)/n-2.0*π,n)] + return [Point{3, Float64}(r*cos(t),r*sin(t),0) for t in range(0.0,(2.0*π)/n-2.0*π,n)] else throw(ArgumentError("Invalid dir specified :$dir, use :acw or :cw")) end end -function circlepoints(f::FunctionType,n; dir=:acw) where {FunctionType <: Function} +function circlepoints(f::FunctionType,n::Int64; dir=:acw) where {FunctionType <: Function} if dir==:acw - return [GeometryBasics.Point{3, Float64}(f(t)*cos(t),f(t)*sin(t),0) for t in range(0,2*π-(2*π)/n,n)] + return [Point{3, Float64}(f(t)*cos(t),f(t)*sin(t),0) for t in range(0,2*π-(2*π)/n,n)] elseif dir==:cw - return [GeometryBasics.Point{3, Float64}(f(t)*cos(t),f(t)*sin(t),0) for t in range(0,(2*π)/n-2*π,n)] + return [Point{3, Float64}(f(t)*cos(t),f(t)*sin(t),0) for t in range(0,(2*π)/n-2*π,n)] end end @@ -1977,7 +1987,7 @@ latter, triangles are formed by slashing the quads. - `V2::Vector`: n-vector """ -function loftlinear(V1,V2;num_steps=nothing,close_loop=true,face_type=:quad) +function loftlinear(V1::Vector{Point{ND,TV}},V2::Vector{Point{ND,TV}};num_steps=nothing,close_loop=true,face_type=:quad) where ND where TV<:Real # Derive num_steps from distance and mean curve point spacing if missing if isnothing(num_steps) @@ -1997,7 +2007,7 @@ function loftlinear(V1,V2;num_steps=nothing,close_loop=true,face_type=:quad) end -function loftpoints2surf(V,num_steps; close_loop=true,face_type=:quad) +function loftpoints2surf(V::Vector{Point{ND,TV}},num_steps; close_loop=true,face_type=:quad) where ND where TV<:Real # Get number of points in each offset curve nc = length(V)/num_steps # Number of points in curve @@ -2096,50 +2106,45 @@ function loftpoints2surf(V,num_steps; close_loop=true,face_type=:quad) end -function dirplot(ax,V,U; color=:black,linewidth=3,scaleval=1.0,style=:from) - E = [GeometryBasics.LineFace{Int}(i,i+length(V)) for i in 1:length(V)] - if style==:from - P = vcat(V,V.+(scaleval.*U)) +function dirplot(ax,V::Vector{Point{ND,TV1}},U::Union{Vector{Point{ND,TV2}},Vector{Vec{ND,TV2}}}; color=:black,linewidth=3,scaleval=1.0,style=:from) where ND where TV1 <: Real where TV2 <: Real + E = [LineFace{Int64}(i,i+length(V)) for i in 1:length(V)] + + if style==:from + P = deepcopy(V) + append!(P,V.+(scaleval.*U)) elseif style==:to - P = vcat(V.-(scaleval.*U),V) + P = V.-(scaleval.*U) + append!(P,V) elseif style==:through UU = (scaleval.*U)/2 - P = vcat(V.-UU,V.+UU) + P = V.-UU + append!(P,V.+UU) else throw(ArgumentError("Invalid style specified :$style, use :from, :to, or :through")) - end + end hp = wireframe!(ax,GeometryBasics.Mesh(P,E),linewidth=linewidth, transparency=false, color=color) return hp end -function normalplot(ax,F::Union{Array{NgonFace{M, Int64}, 1},Array{NgonFace{M, OffsetInteger{-1, UInt32}},1}},V::Vector{Point3{Float64}}; type_flag=:face, color=:black,linewidth=3,scaleval=nothing) where M +function normalplot(ax,F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}; type_flag=:face, color=:black,linewidth=3,scaleval=nothing) where N where TF<:Integer where ND where TV<:Real if isnothing(scaleval) scaleval = pointspacingmean(F,V)/2.0 end if type_flag == :face - N = facenormal(F,V) + NF = facenormal(F,V) V = simplexcenter(F,V) elseif type_flag == :vertex - N = vertexnormal(F,V) + NF = vertexnormal(F,V) else throw(ArgumentError("Incorrect type_flag, use :face or :vertex")) end - hp = dirplot(ax,V,N; color=color,linewidth=linewidth,scaleval=scaleval,style=:from) + hp = dirplot(ax,V,NF; color=color,linewidth=linewidth,scaleval=scaleval,style=:from) return hp end function normalplot(ax,M::GeometryBasics.Mesh; type_flag=:face, color=:black,linewidth=3,scaleval=nothing) F = faces(M) - V = coordinates(M) - - if !isa(F,Array{NgonFace{length(F[1]), Int64}, 1}) - F = togeometrybasics_faces(F) - end - - if !isa(V,Vector{GeometryBasics.Point3{Float64}}) - V = togeometrybasics_points(V) - end - + V = coordinates(M) return normalplot(ax,F,V;type_flag=type_flag, color=color,linewidth=linewidth,scaleval=scaleval) end @@ -2151,7 +2156,7 @@ function wrapindex(i::Int64,n) return 1+mod(i+(n-1),n) end -function edgeangles(F,V) +function edgeangles(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}) where N where TF<:Integer where ND where TV<:Real m = length(F[1]) A = Vector{GeometryBasics.Vec{m, Float64}}() for f in F @@ -2168,12 +2173,12 @@ function edgeangles(F,V) return A end -function quad2tri(F,V; convert_method = :angle)::Vector{TriangleFace{Int64}} +function quad2tri(F::Vector{QuadFace{TF}},V::Vector{Point{ND,TV}}; convert_method = :angle) where TF<:Integer where ND where TV<:Real # Local functions for slash based conversion - forw_slash(f) = [[f[1],f[2],f[3]],[f[3],f[4],f[1]]] # Forward slash - back_slash(f) = [[f[1],f[2],f[4]],[f[2],f[3],f[4]]] # Back slash + forw_slash(f) = [TriangleFace{TF}(f[1],f[2],f[3]),TriangleFace{TF}(f[3],f[4],f[1])] # Forward slash + back_slash(f) = [TriangleFace{TF}(f[1],f[2],f[4]),TriangleFace{TF}(f[2],f[3],f[4])] # Back slash - Ft = Vector{TriangleFace{Int64}}()#(undef,length(Fn1)*2) + Ft = Vector{TriangleFace{TF}}()#(undef,length(Fn1)*2) for f in F if convert_method == :forward ft = forw_slash(f) @@ -2202,10 +2207,10 @@ function quad2tri(F,V; convert_method = :angle)::Vector{TriangleFace{Int64}} return Ft end -function remove_unused_vertices(F,V)::Tuple +function remove_unused_vertices(F,V::Vector{Point{ND,TV}})::Tuple where ND where TV<:Real if isempty(F) # If the face set is empty, return all emtpy outputs Fc = F - Vc = Vector{eltype(V)}(undef,0) + Vc = Vector{Point{ND,TV}}(undef,0) indFix = Vector{Int64}(undef,0) else # Faces not empty to check which indices are used and shorten V if needed indUsed = elements2indices(F) # Indices used @@ -2217,7 +2222,7 @@ function remove_unused_vertices(F,V)::Tuple return Fc, Vc, indFix end -function trisurfslice(F,V,n = Vec{3, Float64}(0.0,1.0,1.0), p = mean(V,dims=1); snapTolerance = 0.0, output_type=:full) +function trisurfslice(F::Vector{TriangleFace{TF}},V::Vector{Point{ND,TV}},n = Vec{3, Float64}(0.0,0.0,1.0), p = mean(V,dims=1); snapTolerance = 0.0, output_type=:full) where TF<:Integer where ND where TV<:Real if !in(output_type,(:full,:above,:below)) throw(ArgumentError("Invalid output_type :$output_type provided, use :full,:above, or :below")) @@ -2231,7 +2236,7 @@ function trisurfslice(F,V,n = Vec{3, Float64}(0.0,1.0,1.0), p = mean(V,dims=1); end LV = d.<0.0 - Fn = Vector{TriangleFace{Int64}}() + Fn = Vector{TriangleFace{TF}}() Cn = Vector{Int64}() Vn = deepcopy(V) D = Dict{Vector{Int64},Int64}() # For pointing from edge to intersection point index @@ -2262,14 +2267,14 @@ function trisurfslice(F,V,n = Vec{3, Float64}(0.0,1.0,1.0), p = mean(V,dims=1); end if output_type == :above || output_type == :full - push!(Fn,TriangleFace{Int64}(D[e1],indP[2],indP[3])) - push!(Fn,TriangleFace{Int64}(D[e1],indP[3],D[e2])) + push!(Fn,TriangleFace{TF}(D[e1],indP[2],indP[3])) + push!(Fn,TriangleFace{TF}(D[e1],indP[3],D[e2])) push!(Cn,1) push!(Cn,1) end if output_type == :below || output_type == :full - push!(Fn,TriangleFace{Int64}(indP[1],D[e1],D[e2])) + push!(Fn,TriangleFace{TF}(indP[1],D[e1],D[e2])) push!(Cn,-1) end @@ -2289,14 +2294,14 @@ function trisurfslice(F,V,n = Vec{3, Float64}(0.0,1.0,1.0), p = mean(V,dims=1); end if output_type == :below || output_type == :full - push!(Fn,TriangleFace{Int64}(D[e1],indP[2],indP[3])) - push!(Fn,TriangleFace{Int64}(D[e1],indP[3],D[e2])) + push!(Fn,TriangleFace{TF}(D[e1],indP[2],indP[3])) + push!(Fn,TriangleFace{TF}(D[e1],indP[3],D[e2])) push!(Cn,-1) push!(Cn,-1) end if output_type == :above || output_type == :full - push!(Fn,TriangleFace{Int64}(indP[1],D[e1],D[e2])) + push!(Fn,TriangleFace{TF}(indP[1],D[e1],D[e2])) push!(Cn,1) end end @@ -2328,14 +2333,14 @@ function count_edge_face(F,E_uni=nothing,indReverse=nothing)::Vector{Int64} return C end -function boundaryedges(F) +function boundaryedges(F::Vector{NgonFace{N,TF}}) where N where TF <: Integer E = meshedges(F) Eu,_,indReverse = gunique(E; return_unique=true, return_index=true, return_inverse=true, sort_entries=true) count_E2F = count_edge_face(F,Eu,indReverse) return Eu[isone.(count_E2F)] end -function edges2curve(Eb) +function edges2curve(Eb::Vector{LineFace{TF}}) where TF <: Integer # TO DO: # Handle while loop safety/breaking # Cope with non-ordered meshes (normals not coherent) @@ -2362,14 +2367,14 @@ end pointspacingmean(F::Array{NgonFace{N, Int64}, 1},V::Vector{Point3{Float64}}) where N The `pointspacingmean` function computes the mean spacing between points. The -input can be just the coordinate set `V`, a vector of GeometryBasics.Point3 +input can be just the coordinate set `V`, a vector of Point3 points, or also a set of edges `E` or faces `F`. If only `V` is provided it is assumed that `V` represents an ordered set of "adjacent" points, e.g. as for a curve. If a vector of edges `E` or a vector of faces `F is also provided, then the average edge length is computed. If instead a set of faces `F` is provided then edges are first computed after which the mean edge spacing is return. """ -function pointspacingmean(V::Vector{Point3{Float64}}) +function pointspacingmean(V::Vector{Point{ND,TV}}) where ND where TV <: Real # Equivalent to: mean(norm.(diff(V,dims=1))) p = 0.0 n = length(V) @@ -2379,7 +2384,7 @@ function pointspacingmean(V::Vector{Point3{Float64}}) return p end -function pointspacingmean(F::Array{NgonFace{N, Int64}, 1},V::Vector{Point3{Float64}}) where N +function pointspacingmean(F::Vector{NgonFace{N,TF}},V::Vector{Point{ND,TV}}) where N where TF<:Integer where ND where TV<:Real if isa(F,Vector{LineFace{Int64}}) E = F else @@ -2398,7 +2403,7 @@ function pointspacingmean(M::GeometryBasics.Mesh) end -function extrudecurve(V1,d; s=1, n=Vec{3, Float64}(0.0,0.0,1.0),num_steps=nothing,close_loop=false,face_type=:quad) +function extrudecurve(V1::Vector{Point{ND,TV}},d; s=1, n=Vec{3, Float64}(0.0,0.0,1.0),num_steps=nothing,close_loop=false,face_type=:quad) where ND where TV<:Real # Derive num_steps from curve point spacing if missing if isnothing(num_steps) num_steps = ceil(Int64,d/pointspacingmean(V1)) @@ -2478,7 +2483,7 @@ function meshgroup(F; con_type = :v) return C end -function distmarch(F,V,indStart; d=nothing, dd=nothing, dist_tol=1e-3,con_V2V=nothing,l=nothing) +function distmarch(F,V::Vector{Point{ND,TV}},indStart; d=nothing, dd=nothing, dist_tol=1e-3,con_V2V=nothing,l=nothing) where ND where TV<:Real # Get vertex-vertex connectivity if isnothing(con_V2V) @@ -2566,12 +2571,12 @@ end # Description This function can compute triangle-ray or triangle-line intersections through -the use of the "Möller-Trumbore triangle-ray intersection algorithm". The +the use of the "Möller-Trumbore triangle-ray intersection algorithm" [1]. The required inputs are as follows: `F` an single face or a vector of faces, e.g. `Vector{TriangleFace{Int64}}` -`V` The triangle vertices as a vector of points, i.e. `Vector{GeometryBasics.Point{3, Float64}}` -`ray_vector` The ray vector which can be `Vector{GeometryBasics.Point{3, Float64}}` or `Vec3{Float64}` +`V` The triangle vertices as a vector of points, i.e. `Vector{Point{3, Float64}}` +`ray_vector` The ray vector which can be `Vector{Point{3, Float64}}` or `Vec3{Float64}` The following optional input parameters can be provided: `rayType = :ray` (default) or `:line`. This defines wether the vector is treated as a ray (extends indefinately) or as a line (finite length) @@ -2582,10 +2587,10 @@ When `triSide=0` both inward and outward intersections are considered. `tolEps = eps(Float64)` (default) # References -[Möller, Tomas; Trumbore, Ben (1997). "Fast, Minimum Storage Ray-Triangle Intersection". Journal of Graphics Tools. 2: 21-28. doi: 10.1080/10867651.1997.10487468.](https://doi.org/10.1080/10867651.1997.10487468) +1. [Möller, Tomas; Trumbore, Ben (1997). "Fast, Minimum Storage Ray-Triangle Intersection". Journal of Graphics Tools. 2: 21-28. doi: 10.1080/10867651.1997.10487468.](https://doi.org/10.1080/10867651.1997.10487468) """ -function ray_triangle_intersect(F::Vector{TriangleFace{Int64}},V,ray_origin,ray_vector; rayType = :ray, triSide = 1, tolEps = eps(Float64)) - P = Vector{GeometryBasics.Point{3, Float64}}() +function ray_triangle_intersect(F::Vector{TriangleFace{TF}},V::Vector{Point{ND,TV}},ray_origin,ray_vector; rayType = :ray, triSide = 1, tolEps = eps(Float64)) where TF <: Integer where ND where TV<:Real + P = Vector{Point{ND,TV}}() indIntersect = Vector{Int64}() for qf in eachindex(F) p = ray_triangle_intersect(F[qf],V,ray_origin,ray_vector; rayType = rayType, triSide = triSide, tolEps = tolEps) @@ -2597,7 +2602,7 @@ function ray_triangle_intersect(F::Vector{TriangleFace{Int64}},V,ray_origin,ray_ return P,indIntersect end -function ray_triangle_intersect(f::TriangleFace{Int64},V,ray_origin,ray_vector; rayType = :ray, triSide = 1, tolEps = eps(Float64)) +function ray_triangle_intersect(f::TriangleFace{Int64},V::Vector{Point{ND,TV}},ray_origin,ray_vector; rayType = :ray, triSide = 1, tolEps = eps(Float64)) where ND where TV<:Real # Edge vectors P1 = V[f[1]] # First corner point @@ -2615,7 +2620,7 @@ function ray_triangle_intersect(f::TriangleFace{Int64},V,ray_origin,ray_vector; boolDet = det_vec reverse(f),F) # [NgonFace{N, Int64}(reverse(f)) for f in F] end """ - R = kabsch_rot(V1::Array{Point{N, Float64}, 1},V2::Array{Point{N, Float64}, 1}) where N + R = kabsch_rot(V1::Array{Point{N, T}, 1},V2::Array{Point{N, TT}, 1}) where N where T<:Real where TT<:Real # Description Computes the rotation tensor `R` to rotate the points in `V1` to best match the @@ -2818,7 +2823,7 @@ points in `V2`. [Wolfgang Kabsch, A solution for the best rotation to relate two sets of vectors, Acta Crystallographica Section A, vol. 32, no. 5, pp. 922-923, 1976, doi: 10.1107/S0567739476001873](https://doi.org/10.1107/S0567739476001873) [https://en.wikipedia.org/wiki/Kabsch_algorithm](https://en.wikipedia.org/wiki/Kabsch_algorithm) """ -function kabsch_rot(V1::Array{Point{N, Float64}, 1},V2::Array{Point{N, Float64}, 1}) where N +function kabsch_rot(V1::Vector{Point{ND,TV1}},V2::Vector{Point{ND,TV2}}) where ND where TV1<:Real where TV2<:Real # Centre on means V1 = V1.-mean(V1) V2 = V2.-mean(V2) @@ -2857,7 +2862,7 @@ number (negative or positive) of full twists to the loft. Finally the optional parameter `close_loop` (default is `true`) determines if the section curves are deemed closed or open ended. """ -function sweeploft(Vc,V1,V2; face_type=:quad, num_twist = 0, close_loop=true) +function sweeploft(Vc::Vector{Point{ND,TV}},V1::Vector{Point{ND,TV}},V2::Vector{Point{ND,TV}}; face_type=:quad, num_twist = 0, close_loop=true) where ND where TV<:Real np = length(V1) # Number of section points nc = length(Vc) # Number of curve steps @@ -2907,7 +2912,7 @@ function sweeploft(Vc,V1,V2; face_type=:quad, num_twist = 0, close_loop=true) end if q == nc-1 # Once here, second to last, a potential rotational mismatch needs to be resolve - Q_fix = Rotations.RotMatrix3{Float64}(S2\S2p) # Rotation between last and second to last + Q_fix = RotMatrix3{Float64}(S2\S2p) # Rotation between last and second to last t_a = Rotations.params(AngleAxis(Q_fix)) # Angle/axis representation if dot(t_a[2:end],n3)<0 β_fix = t_a[1] @@ -2933,7 +2938,7 @@ function sweeploft(Vc,V1,V2; face_type=:quad, num_twist = 0, close_loop=true) end -function revolvecurve(Vc,θ=2.0*pi; s=0, n=Vec{3, Float64}(0.0,0.0,1.0),num_steps=nothing,close_loop=false,face_type=:quad) +function revolvecurve(Vc::Vector{Point{ND,TV}},θ=2.0*pi; s=0, n=Vec{3, Float64}(0.0,0.0,1.0),num_steps=nothing,close_loop=false,face_type=:quad) where ND where TV<:Real # Compute num_steps from curve point spacing if isnothing(num_steps) rMax = 0.0 @@ -2967,4 +2972,64 @@ function revolvecurve(Vc,θ=2.0*pi; s=0, n=Vec{3, Float64}(0.0,0.0,1.0),num_step end return loftpoints2surf(V,num_steps;close_loop=close_loop,face_type=face_type) +end + + +""" + batman(n::Int64) +# Description +The `batman` function creates points on the curve for the Batman logo. The curve +is useful for testing surface meshing algorithms since it contains sharp +transitions and pointy features. The user requests `n` points on the curve. The +default forces exactly `n` points which may result in an assymetric curve. To +instead force symmetry the user can set the optional parameter `symmetric=true`. +In this case the output will be symmetric allong the y-axis, however the number +of points on the curve may have increased (if the input `n` is not even). The +second optional input is the direction of the curve, i.e. if it is clockwise, +`dir=:cw` or anti-clockwise `dir=:acw`. +The implementation is based on a "parameterised Batman equation" [1](https://www.desmos.com/calculator/ajnzwedvql). +The following modifications where made, the curve is here centered around +[0,0,0], scaled to be 2 in width, resampled evenly, and the default curve +direction is anti-clockwise. + +# References +1. https://www.desmos.com/calculator/ajnzwedvql +""" +function batman(n::Int64; symmetric = false, dir=:acw) + tt = range(8,0.0,n) + + x = [( 0.3*t + 0.2*abs(t-1.0) + 2.2*abs(t-2.0) - 2.7*abs(t-3.0) -3.0*abs(t-5.0) + 3*abs(t-7.0) + + 5.0*sin(π/4.0*(abs(t-3.0)-abs(t-4.0)+1.0)) + + 5.0/4.0 * (abs(t-4.0)-abs(t-5.0)-1.0)^3.0 + - 5.3*cos((π/2.0 + asin(47/53))*(abs(t-7.0) + - abs(t-8.0) - 1.0)/2.0) + 2.8 ) for t in tt] + + y = [( 3.0/2.0*abs(t-1.0) - 3.0/2.0*abs(t-2.0) - 29.0/4.0*abs(t-4.0) + 29.0/4.0*abs(t-5.0) + + 7.0/16.0*(abs(t-2.0) - abs(t-3.0) -1.0)^4.0 + + 4.5*sin( π/4.0*(abs(t-3.0) - abs(t-4.0) - 1.0) ) + - ((3.0/5.0*sqrt(2)) * abs(abs(t-5.0)-abs(t-7.0))^(5.0/2.0)) + + 6.4 * sin( ( (π/2.0 + asin(47/53)) * ((abs(t-7.0) + - abs(t-8.0) + 1.0)/2.0) ) + asin(56/64) ) + 4.95 ) for t in tt] + + if symmetric == true + if iseven(n) + m = round(Int64,n/2+1) + else + m = ceil(Int64,n/2)+1 + end + V,_ = evenly_sample([Point{3, Float64}(x[i]/22,y[i]/22,0.0) for i in eachindex(x)],m) + V2 = [Point{3, Float64}(-V[i][1],V[i][2],0.0) for i in length(V)-1:-1:2] + append!(V,V2) + else + V = [Point{3, Float64}(x[i]/22,y[i]/22,0.0) for i in eachindex(x)] + V2 = [Point{3, Float64}(-V[i][1],V[i][2],0.0) for i in length(V)-1:-1:2] + append!(V,V2) + V,_ = evenly_sample(V,n) + end + if dir==:cw + reverse!(V) + circshift!(V,1) + end + V = V.-mean(V,dims=1) + return V end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 045d253..27b001b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -336,24 +336,6 @@ end @test isapprox(V, expected, atol = eps_level) end - @testset "Vector Vector" begin - P = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [1.0, 1.0, 0.0], [1.0, 1.0, 1.0]] - n = 10 # Number of points - V = nbezier(P,n) # Get Bezier fit points - expected = [[0.0, 0.0, 0.0], - [0.29766803840877915, 0.034293552812071325, 0.0013717421124828531], - [0.5294924554183813, 0.1262002743484225, 0.010973936899862825], - [0.7037037037037037, 0.25925925925925924, 0.037037037037037035], - [0.8285322359396433, 0.4170096021947874, 0.0877914951989026], - [0.9122085048010974, 0.5829903978052127, 0.1714677640603567], - [0.9629629629629629, 0.7407407407407407, 0.2962962962962963], - [0.9890260631001372, 0.8737997256515775, 0.4705075445816187], - [0.9986282578875172, 0.9657064471879286, 0.7023319615912208], - [1.0, 1.0, 1.0]] - @test typeof(V) == typeof(P) - @test length(V) == n - @test isapprox(V, expected, atol = eps_level) - end end @@ -870,7 +852,11 @@ end end -@testset "togeometrybasics_faces" verbose = true begin +@testset "tofaces" verbose = true begin + # Edges matrix and vector + Fem = [1 2; 4 5] + Fev = [[1,2],[4,5]] + # Triangles matrix and vector Ftm = [1 2 3; 4 5 6] Ftv = [[1,2,3],[4,5,6]] @@ -887,84 +873,96 @@ end fileName_mesh = joinpath(comododir(),"assets","stl","stanford_bunny_low.stl") M = load(fileName_mesh) - @testset "Matrix input" begin - Ftm_G = togeometrybasics_faces([Ftm[1,:]]) + @testset "Matrix input" verbose = true begin + Fem_G = tofaces([Fem[1,:]]) + @testset "1 LineFace" begin + @test isa(Fem_G,Vector{GeometryBasics.LineFace{Int64}}) + @test length(Fem_G) == 1 + end + + Ftm_G = tofaces([Ftm[1,:]]) @testset "1 TriangleFace" begin @test isa(Ftm_G,Vector{GeometryBasics.TriangleFace{Int64}}) @test length(Ftm_G) == 1 end - Ftm_G = togeometrybasics_faces(Ftm) + Ftm_G = tofaces(Ftm) @testset "TriangleFace" begin @test isa(Ftm_G,Vector{GeometryBasics.TriangleFace{Int64}}) @test length(Ftm_G) == size(Ftm,1) end - Fqm_G = togeometrybasics_faces(Fqm) + Fqm_G = tofaces(Fqm) @testset "QuadFace" begin @test isa(Fqm_G,Vector{GeometryBasics.QuadFace{Int64}}) @test length(Fqm_G) == size(Fqm,1) end - Fnm_G = togeometrybasics_faces(Fnm) + Fnm_G = tofaces(Fnm) @testset "NgonFace" begin @test isa(Fnm_G,Vector{GeometryBasics.NgonFace{5,Int64}}) @test length(Fnm_G) == size(Fnm,1) end end - @testset "vector input" begin - Ftv_G = togeometrybasics_faces([Ftv[1]]) + @testset "vector input" verbose = true begin + Fev_G = tofaces([Fev[1]]) + @testset "1 LineFace" begin + @test isa(Fev_G,Vector{GeometryBasics.LineFace{Int64}}) + @test length(Fev_G) == 1 + end + + Ftv_G = tofaces([Ftv[1]]) @testset "1 TriangleFace" begin @test isa(Ftv_G,Vector{GeometryBasics.TriangleFace{Int64}}) @test length(Ftv_G) == 1 end - Ftv_G = togeometrybasics_faces(Ftv) + Ftv_G = tofaces(Ftv) @testset "TriangleFace" begin @test isa(Ftv_G,Vector{GeometryBasics.TriangleFace{Int64}}) @test length(Ftv_G) == length(Ftv) end - Fqv_G = togeometrybasics_faces(Fqv) + Fqv_G = tofaces(Fqv) @testset "QuadFace" begin @test isa(Fqv_G,Vector{GeometryBasics.QuadFace{Int64}}) @test length(Fqv_G) == length(Fqv) end - Fnv_G = togeometrybasics_faces(Fnv) + Fnv_G = tofaces(Fnv) @testset "NgonFace" begin @test isa(Fnv_G,Vector{GeometryBasics.NgonFace{5,Int64}}) @test length(Fnv_G) == length(Fnv) - end - - F = togeometrybasics_faces(faces(M)) - @testset "Vector{NgonFace{3, OffsetInteger{-1, UInt32}}}" begin - @test isa(F,Vector{GeometryBasics.TriangleFace{Int64}}) - @test length(F) == length(faces(M)) - end + end + + Fnv_G2 = tofaces(Fnv_G) + @testset "NgonFace no change" begin + @test isa(Fnv_G2,Vector{GeometryBasics.NgonFace{5,Int64}}) + @test length(Fnv_G2) == length(Fnv_G) + end end end -@testset "togeometrybasics_points" verbose = true begin +@testset "topoints" verbose = true begin @testset "Matrix input" begin - V = togeometrybasics_points(rand(10,3)) + V = topoints(rand(10,3)) @test isa(V,Vector{GeometryBasics.Point3{Float64}}) @test length(V) == 10 end @testset "Vector Float64" begin Vv = [rand(3) for _ in 1:5] - V = togeometrybasics_points(Vv) + V = topoints(Vv) @test isa(V,Vector{GeometryBasics.Point3{Float64}}) @test length(V) == 5 end @testset "Vector Vec3" begin Vv = Vector{Vec3{Float64}}(undef,5) - V = togeometrybasics_points(Vv) + V = topoints(Vv) @test isa(V,Vector{GeometryBasics.Point3{Float64}}) @test length(V) == 5 end @@ -972,13 +970,13 @@ end @testset "Vector Vec{m,Float64}}" begin m = 4 Vv = Vector{Vec{m,Float64}}(undef,5) - V = togeometrybasics_points(Vv) + V = topoints(Vv) @test isa(V,Vector{GeometryBasics.Point{m,Float64}}) @test length(V) == 5 m = 5 Vv = Vector{Vec{m,Float64}}(undef,5) - V = togeometrybasics_points(Vv) + V = topoints(Vv) @test isa(V,Vector{GeometryBasics.Point{m,Float64}}) @test length(V) == 5 end @@ -988,12 +986,22 @@ end fileName_mesh = joinpath(comododir(),"assets","stl","stanford_bunny_low.stl") M = load(fileName_mesh) Vv = coordinates(M) - V = togeometrybasics_points(Vv) + V = topoints(Vv) + @test isa(V,Vector{GeometryBasics.Point3{Float32}}) + @test length(V) == length(Vv) + end + + @testset "points no change" begin + # Imported triangular mesh + Vv = rand(Point{3,Float64},5) + V = topoints(Vv) @test isa(V,Vector{GeometryBasics.Point3{Float64}}) @test length(V) == length(Vv) end + end + @testset "togeometrybasics_mesh" verbose = true begin @testset "Matrix input" begin @@ -1256,17 +1264,6 @@ end @test edgelengths(F,V*pi) == pi.*[1.0, 1.0, 1.0, 1.0] # Scaled square end - @testset "F::Vector{Vector{Int64}}, V::Vector{Vec3}" begin - F = [[1,2,3,4]] - V = Vector{Vec3{Float64}}(undef,4) - V[1] = Vec3(0.0, 0.0, 0.0) - V[2] = Vec3(1.0, 0.0, 0.0) - V[3] = Vec3(1.0, 1.0, 0.0) - V[4] = Vec3(0.0, 1.0, 0.0) - - @test edgelengths(F,V) == [1.0, 1.0, 1.0, 1.0] - end - @testset "GeometryBasics LineFace edges" begin F = [QuadFace{Int64}(1, 2, 3, 4)] E = meshedges(F) @@ -2448,6 +2445,7 @@ end end end + @testset "normalizevector" verbose = true begin eps_level = 1e-4 @@ -2457,7 +2455,7 @@ end @test n isa typeof(v) @test n == [0.0, 0.0, 1.0] - v = [1.0, 1.0, 0.0] + v = Vec{3,Float64}(1.0, 1.0, 0.0) n = normalizevector(v) @test n isa typeof(v) @test isapprox(n,√2/2*[1.0, 1.0, 0.0], atol=eps_level) @@ -2473,17 +2471,15 @@ end F = faces(M) V = coordinates(M) N = facenormal(F,V) - U = [n./(1.0.+rand(1)) for n in N] + U = [Vec{3,Float64}(n./(1.0.+rand(1))) for n in N] NN = normalizevector(U) @test eltype(NN) == eltype(U) @test isapprox(NN,N, atol=eps_level) end - end - @testset "circlepoints" verbose = true begin eps_level = 1e-4 @@ -2794,7 +2790,6 @@ end @test typeof(hp1) == Wireframe{Tuple{GeometryBasics.Mesh{3, Float64, Line{3, Float64}, SimpleFaceView{3, Float64, 2, Int64, Point3{Float64}, LineFace{Int64}}}}} @test length(faces(Mp)) == length(V) - hp1 = normalplot(ax,F,V; type_flag=:vertex, color=:black,linewidth=3,scaleval=nothing) Mp = hp1[1].val @test typeof(hp1) == Wireframe{Tuple{GeometryBasics.Mesh{3, Float64, Line{3, Float64}, SimpleFaceView{3, Float64, 2, Int64, Point3{Float64}, LineFace{Int64}}}}} @@ -2803,10 +2798,11 @@ end fileName_mesh = joinpath(comododir(),"assets","obj","spot_control_mesh.obj") M = load(fileName_mesh) F = faces(M) - V = coordinates(M) - hp1 = normalplot(ax,M; type_flag=:vertex, color=:black,linewidth=3,scaleval=nothing) + V = topoints(coordinates(M)) + + hp1 = normalplot(ax,F,V; type_flag=:vertex, color=:black,linewidth=3,scaleval=nothing) Mp = hp1[1].val - @test typeof(hp1) == Wireframe{Tuple{GeometryBasics.Mesh{3, Float64, Line{3, Float64}, SimpleFaceView{3, Float64, 2, Int64, Point3{Float64}, LineFace{Int64}}}}} + @test typeof(hp1) == Wireframe{Tuple{GeometryBasics.Mesh{3, Float32, Line{3, Float32}, SimpleFaceView{3, Float32, 2, Int64, Point3{Float32}, LineFace{Int64}}}}} @test length(faces(Mp)) == length(V) # Not supported yet @@ -2875,7 +2871,7 @@ end fDef[1,2] = tan(a) # Sheared cube coordinates - V2 = togeometrybasics_points([fDef*v for v in V]) + V2 = topoints([fDef*v for v in V]) A = edgeangles(F,V) # Angles for regular cube A2 = edgeangles(F,V2) # Angles for sheared cube @@ -2896,7 +2892,7 @@ end end a = pi/6 # "45 degree shear" f[1,2] = tan(a) - V = togeometrybasics_points([f*v for v in V]) # Shear the cube + V = topoints([f*v for v in V]) # Shear the cube # Build deformation gradient tensor to induce shear with known angles f = zeros(3,3) @@ -2905,7 +2901,7 @@ end end a = pi/6 # "45 degree shear" f[3,1] = tan(a) - V = togeometrybasics_points([f*v for v in V]) # Shear the cube + V = topoints([f*v for v in V]) # Shear the cube @testset "Errors" begin @test_throws ArgumentError quad2tri(F,V; convert_method = :wrong) @@ -3148,19 +3144,19 @@ end end @testset "Direction (s) variations" begin - 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, d; s=1, n=Vec{3, Float64}(0.0,0.0,1.0), num_steps=num_steps, close_loop=true, face_type=:quad) z = [v[3] for v in V] zMax = maximum(z) zMin = minimum(z) @test isapprox(zMax,d,atol = eps_level) && isapprox(zMin,0.0,atol = eps_level) - F, V = extrudecurve(Vc, d; s=0, n=[0.0,0.0,1.0], num_steps=num_steps, close_loop=true, face_type=:quad) + F, V = extrudecurve(Vc, d; s=0, n=Vec{3, Float64}(0.0,0.0,1.0), num_steps=num_steps, close_loop=true, face_type=:quad) z = [v[3] for v in V] zMax = maximum(z) zMin = minimum(z) @test isapprox(zMax,d/2,atol = eps_level) && isapprox(zMin,-d/2,atol = eps_level) - 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, d; s=-1, n=Vec{3, Float64}(0.0,0.0,1.0), num_steps=num_steps, close_loop=true, face_type=:quad) z = [v[3] for v in V] zMax = maximum(z) zMin = minimum(z) @@ -3168,21 +3164,21 @@ end end @testset "Direction (n) variations" begin - n=[0.0,0.0,1.0] # Upward + n=Vec{3, Float64}(0.0,0.0,1.0) # Upward F, V = extrudecurve(Vc, d; s=1, n=n, num_steps=num_steps, close_loop=true, face_type=:quad) z = [v[3] for v in V] zMax = maximum(z) zMin = minimum(z) @test isapprox(zMax,d,atol = eps_level) && isapprox(zMin,0.0,atol = eps_level) - n=[0.0,0.0,-1.0] # Downward + n=Vec{3, Float64}(0.0,0.0,-1.0) # Downward F, V = extrudecurve(Vc, d; s=1, n=n, num_steps=num_steps, close_loop=true, face_type=:quad) z = [v[3] for v in V] zMax = maximum(z) zMin = minimum(z) @test isapprox(zMax,0.0,atol = eps_level) && isapprox(zMin,-d,atol = eps_level) - n = normalizevector([1.0,0.0,1.0]) # 45 degree direction upward + n = normalizevector(Vec{3, Float64}(1.0,0.0,1.0)) # 45 degree direction upward F, V = extrudecurve(Vc, d; s=1, n=n, num_steps=num_steps, close_loop=true, face_type=:quad) z = [v[3] for v in V] zMax = maximum(z) @@ -3191,7 +3187,7 @@ end end @testset "face_type=:quad" begin - 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, d; s=1, n=Vec{3, Float64}(0.0,0.0,1.0), num_steps=num_steps, close_loop=true, face_type=:quad) z = [v[3] for v in V] zMax = maximum(z) zMin = minimum(z) @@ -3208,7 +3204,7 @@ end end @testset "face_type=:tri" begin - F, V = extrudecurve(Vc, d; s=1, n=[0.0,0.0,1.0], num_steps=num_steps, close_loop=true, face_type=:tri) + F, V = extrudecurve(Vc, d; s=1, n=Vec{3, Float64}(0.0,0.0,1.0), num_steps=num_steps, close_loop=true, face_type=:tri) z = [v[3] for v in V] zMax = maximum(z) zMin = minimum(z) @@ -3225,7 +3221,7 @@ end end @testset "face_type=:tri_slash" begin - F, V = extrudecurve(Vc, d; s=1, n=[0.0,0.0,1.0], num_steps=num_steps, close_loop=true, face_type=:tri_slash) + F, V = extrudecurve(Vc, d; s=1, n=Vec{3, Float64}(0.0,0.0,1.0), num_steps=num_steps, close_loop=true, face_type=:tri_slash) z = [v[3] for v in V] zMax = maximum(z) zMin = minimum(z) @@ -3242,7 +3238,7 @@ end end @testset "face_type=:quad2tri" begin - F, V = extrudecurve(Vc, d; s=1, n=[0.0,0.0,1.0], num_steps=num_steps, close_loop=true, face_type=:quad2tri) + F, V = extrudecurve(Vc, d; s=1, n=Vec{3, Float64}(0.0,0.0,1.0), num_steps=num_steps, close_loop=true, face_type=:quad2tri) z = [v[3] for v in V] zMax = maximum(z) zMin = minimum(z) @@ -4024,3 +4020,46 @@ end end end + +@testset "batman" verbose=true begin + eps_level = 1e-6 + + n = 50 # Number of points, first testig an even number + V = batman(n) + ind = round.(Int64,range(1,length(V),5)) + + @test length(V) == n # Correct length + @test isa(V,Vector{Point{3,Float64}}) # Correct type + # Correct coordinates for this case + @test isapprox(V[ind],Point3{Float64}[ + [4.959248144765713e-5, -0.6888052278988526, 0.0], + [0.9989308524502467, -0.12980802206407838, 0.0], + [-0.004310373406159424, 0.3627431045478119, 0.0], + [-0.9715858301799283, -0.014868116101481532, 0.0], + [-0.022332909569322653, -0.5717894489374953, 0.0]],atol=eps_level) + + V = batman(n; dir=:cw) + @test length(V) == n # Correct length + @test isa(V,Vector{Point{3,Float64}}) + + V = batman(n; dir=:acw, symmetric=false) + @test length(V) == n # Correct length + @test isa(V,Vector{Point{3,Float64}}) + + V = batman(n; dir=:acw, symmetric=true) + @test length(V) == n # Correct length + @test isa(V,Vector{Point{3,Float64}}) + + m = n+1 # force uneven + V = batman(m; dir=:acw, symmetric=true) + @test length(V) == m+1 # Correct length + @test isa(V,Vector{Point{3,Float64}}) + + V = batman(m; dir=:cw, symmetric=true) + @test length(V) == m+1 # Correct length + @test isa(V,Vector{Point{3,Float64}}) + + +end + +