Fixed issue in regiontrimesh

COMODO-research · Feb 13, 2025 · 5e368bf · 5e368bf
1 parent b1d5f7f
commit 5e368bf
Show file tree

Hide file tree

Showing 6 changed files with 131 additions and 65 deletions.
diff --git a/Project.toml b/Project.toml
@@ -19,7 +19,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TetGen = "c5d3f3f7-f850-59f6-8a2e-ffc6dc1317ea"
 
 [compat]
-BSplineKit = "0.17"
+BSplineKit = "0.18"
 DataStructures = "0.18.16"
 DelaunayTriangulation = "1.6.3"
 Distances = "0.10"

diff --git a/examples/demo_regiontrimesh.jl b/examples/demo_regiontrimesh.jl
@@ -7,7 +7,7 @@ This demo shows the use of `hexbox` to generate a hexahedral mesh for a 3D box
 domain. 
 =#
 
-testCase = 5
+testCase = 6
 
 if testCase == 1 # Batman curve
     n = 120
@@ -37,7 +37,36 @@ elseif testCase == 3
     VT = (V,)
     R = ([1],)
     P = (pointSpacing)
-elseif testCase == 4
+elseif testCase == 4 
+
+    function squarepoints(w,h,pointSpacing; dir=:acw)
+        nw = ceil(Int,w./pointSpacing)+1
+        nh = ceil(Int,h./pointSpacing)+1
+        sw = w./(nw-1)
+        sh = h./(nh-1)        
+        if dir == :acw            
+            V1 = collect(range(Point{3,Float64}( w./2.0,           h/2.0, 0.0), Point{3,Float64}(-w./2.0,         h/2.0, 0.0),   nw))
+            V2 = collect(range(Point{3,Float64}(-w./2.0,      h/2.0 - sh, 0.0), Point{3,Float64}(-w./2.0,        -h/2.0, 0.0), nh-1))
+            V3 = collect(range(Point{3,Float64}(-w./2.0+sw,       -h/2.0, 0.0), Point{3,Float64}( w./2.0,        -h/2.0, 0.0), nw-1))
+            V4 = collect(range(Point{3,Float64}( w./2.0,     -h/2.0 + sh, 0.0), Point{3,Float64}( w./2.0,     h/2.0 -sh, 0.0), nh-2))            
+        else
+            V1 = collect(range(Point{3,Float64}(-w./2.0,           h/2.0, 0.0), Point{3,Float64}( w./2.0,         h/2.0, 0.0),   nw))
+            V2 = collect(range(Point{3,Float64}( w./2.0,      h/2.0 - sh, 0.0), Point{3,Float64}( w./2.0,        -h/2.0, 0.0), nh-1))
+            V3 = collect(range(Point{3,Float64}( w./2.0-sw,       -h/2.0, 0.0), Point{3,Float64}(-w./2.0,        -h/2.0, 0.0), nw-1))
+            V4 = collect(range(Point{3,Float64}(-w./2.0,     -h/2.0 + sh, 0.0), Point{3,Float64}(-w./2.0,     h/2.0 -sh, 0.0), nh-2))            
+        end
+        return [V1; V2; V3; V4]        
+    end
+
+    w = 4.0
+    h = 8.0
+    pointSpacing = 1.0
+    V = squarepoints(w,h,pointSpacing; dir=:acw)
+
+    VT = (V,)
+    R = ([1],)
+    P = (pointSpacing)
+elseif testCase == 5
     rFun1(t) = 12.0 + 3.0.*sin(5*t)
     n1 = 120
     V1 = circlepoints(rFun1,n1)
@@ -61,18 +90,17 @@ elseif testCase == 4
     VT = (V1,V2,V3,V4,V5,V6) # Curves
     R = ([1,2],[2,3,4,5,6],[4],[5]) # Regions 
     P = (1,0.6,0.2,0.3)  # Point spacings
-
-elseif testCase == 5    
-    n1 = 120
+elseif testCase == 6    
+    n1 = 100
     r1 = 20.0
     V1 = circlepoints(r1,n1)
 
-    n2 = 100
+    n2 = 90
     r2 = 12.0
     V2 = circlepoints(r2,n2)
     V2 = [Point{3,Float64}(v[1]+6,v[2],v[3]) for v in V2]
 
-    n3 = 30
+    n3 = 25
     r3 = 4
     V3 = circlepoints(r3,n3)
     V3 = [Point{3,Float64}(v[1]-14,v[2],v[3]) for v in V3]
@@ -89,11 +117,12 @@ elseif testCase == 5
 
     VT = (V1,V2,V3,V4,V5) # Curves
     R = ([1,2,3],[2,4,5],[5]) # Regions 
-    P = (1,0.75,0.5)  # Point spacings
+    P = (1.5,0.75,0.5)  # Point spacings
 end
 
-F,V,C = regiontrimesh(VT,R,P)
+VTp = deepcopy(VT) # Copy for plotting
 
+F,V,C = regiontrimesh(VT,R,P)
 
 # Visualisation
 Fp,Vp = separate_vertices(F,V) # Give each face its own point set 
@@ -102,5 +131,9 @@ Cp = simplex2vertexdata(Fp,C) # Convert face color data to vertex color data
 fig = Figure(size=(1200,1000))
 ax1 = Axis3(fig[1, 1],aspect = :data,title="Multi-region meshing",azimuth=-pi/2,elevation=pi/2)
 hp1 = poly!(ax1,GeometryBasics.Mesh(Vp,Fp), strokewidth=1,color=Cp, strokecolor=:black, shading = FastShading, transparency=false,colormap=Makie.Categorical(Makie.Reverse(:Spectral)))
+# for V in VTp
+#     scatter!(ax1, V, markersize=15, color=:black)
+#     lines!(ax1, V, linewidth=5, color=:black)
+# end
 Colorbar(fig[1, 1][1, 2], hp1)
 fig
diff --git a/examples/demo_remove_unused_vertices.jl b/examples/demo_remove_unused_vertices.jl
@@ -2,22 +2,40 @@ using Comodo
 using Comodo.GLMakie
 using Comodo.GeometryBasics
 
-# Example geometry
-F,V = geosphere(3,1.0)
-VC = simplexcenter(F,V)
-F = [F[i] for i in findall(map(v-> v[3]>0,VC))] # Remove some faces
+testCase = 1
+if testCase == 1
+    # Example geometry
+    F,V = geosphere(3,1.0)
+    VC = simplexcenter(F,V)
+    F = [F[i] for i in findall(map(v-> v[3]>0,VC))] # Remove some faces
+elseif testCase == 2
+    F,V = geosphere(3,1.0)
+    Vr = 2.0*randn(Point{3,Float64},5)
+    V = [Vr;V]
+    F = [f.+length(Vr) for f in F]
+elseif testCase == 3
+    F,V = geosphere(2,1.0)
+    F = [F[i] for i in 1:2:length(F)]
+end
 
-Fc,Vc = remove_unused_vertices(F,V)
+indCheck = intersect(findall([v[1]>0 for v in V]),reduce(vcat,F))
+Vn = deepcopy(V)
+
+Fc,Vc,indFix = remove_unused_vertices(F,V)
+indCheck_c = indFix[indCheck]
 
 # Visualization
-markersize = 20
-M = GeometryBasics.Mesh(Vc,Fc)
+markersize = 10
+markersize2 = 20
 fig = Figure(size=(1200,1200))
-ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "A cut mesh with unused vertices removed")
+ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Mesh with unused vertices")
+hp1 = poly!(ax1,GeometryBasics.Mesh(V,F), strokewidth=1,color=:white, strokecolor=:blue, shading = FastShading, transparency=false)
+scatter!(Vn,markersize=markersize,color=:red)
+scatter!(Vn[indCheck],markersize=markersize2,color=:blue)
 
-hp3 = poly!(ax1,M, strokewidth=1,color=:white, strokecolor=:blue, shading = FastShading, transparency=false)
-# hp3 = normalplot(ax1,M)
+ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Mesh with unused vertices removed")
+hp2 = poly!(ax2,GeometryBasics.Mesh(Vc,Fc), strokewidth=1,color=:white, strokecolor=:blue, shading = FastShading, transparency=false)
 scatter!(Vc,markersize=markersize,color=:red)
-
+scatter!(Vc[indCheck_c],markersize=markersize2,color=:blue)
 
 fig
diff --git a/src/Comodo.jl b/src/Comodo.jl
@@ -69,7 +69,7 @@ export squircle, circlerange, edgefaceangles, faceanglesegment, eulerchar
 export rhombicdodecahedronfoam, kelvinfoam, truncatedoctahedron, ntrapezohedron, hexagonaltrapezohedron #, tetrakaidecahedron
 export mag, indexmap!, indexmap, minp, maxp, spacing2numvertices
 export joingeom, quadbox, tribox, tetbox, pad3, getisosurface 
-export randangle, stepfunc, perlin_noise
+export randangle, stepfunc, perlin_noise, removepoints
 
 # Export functions: Visualization related
 export slidercontrol, slider2anim, dirplot, normalplot

diff --git a/src/functions.jl b/src/functions.jl
@@ -3131,9 +3131,9 @@ end
 function remove_unused_vertices(F,V::Vector{Point{ND,TV}}) where ND where TV<:Real
     if isempty(F) # If the face set is empty, return all emtpy outputs
         Fc = F
-        Vc = Vector{Point{ND,TV}}(undef,0)
-        indFix = Vector{Int}(undef,0)
-    else # Faces not empty to check which indices are used and shorten V if needed        
+        Vc = Vector{Point{ND,TV}}()
+        indFix = Vector{Int}()
+    else # Faces not empty, so check which indices are used and shorten V if needed        
         indUsed = elements2indices(F) # Indices used
         Vc = V[indUsed] # Remove unused points    
         indFix = zeros(Int,length(V))
@@ -4405,8 +4405,7 @@ function regiontrimesh(VT,R,P)
         r = R[q] # The curve indices for the current region       
         pointSpacing = P[q] # Region point spacing
         Vn = reduce(vcat,VT[r]) # The current curve point set 
-        numCurvePoints = length(Vn) # Number of points in the current curve set
-
+
         # Creating constraints
         constrained_segments = Vector{Vector{Vector{Int}}}()
         n = 1        
@@ -4441,53 +4440,46 @@ function regiontrimesh(VT,R,P)
         constrained_segments_ori = deepcopy(constrained_segments) # Clone since triangulate can add new constraint points
         TRn = triangulate(Vn; boundary_nodes=constrained_segments, delete_ghosts=true)
         Fn = [TriangleFace{Int}(tr) for tr in each_solid_triangle(TRn)] 
-
+        Vn = get_points(TRn)
+        # Fn,Vn,indFix = remove_unused_vertices(Fn,Vn)
+        # constrained_segments = [[indFix[c[1]]] for c in constrained_segments_ori] # Fix indices after point removal 
+
         # Check if new boundary points were introduced and remove if needed 
         Eb = boundaryedges(Fn)
         indB = unique(reduce(vcat,Eb))
-        indRemove = indB[indB.>numCurvePoints]        
-        if !iszero(length(indRemove))                        
-            logicKeep = fill(true,length(Vn))
-            logicKeep[indRemove] .= false
-            indKeep = findall(logicKeep)
-            indFix = zeros(Int,length(Vn))
-            indFix[indKeep] = 1:length(indKeep)
-            Vn = Vn[indKeep]
-            constrained_segments = [[indFix[c[1]]] for c in constrained_segments_ori] # Fix indices after point removal                
+        indConstrained = reduce(vcat,reduce(vcat,constrained_segments_ori))
+        indRemove = setdiff(indB,indConstrained)
+        if !isempty(indRemove)    
+            Vn,indFix = removepoints(Vn,indRemove)
+            constrained_segments = [[indFix[c[1]]] for c in constrained_segments_ori]
+
+            # Redo triangulation after points have been removed
+            constrained_segments_ori = deepcopy(constrained_segments) # Clone since triangulate can add new constraint points
             TRn = triangulate(Vn; boundary_nodes=constrained_segments,delete_ghosts=true)
             Fn = [TriangleFace{Int}(tr) for tr in each_solid_triangle(TRn)] 
             Vn = get_points(TRn)
+            constrained_segments = constrained_segments_ori
         end
-
-        # Remove unused points (e.g. outside region)
-        Fn,Vn,indFix1 = remove_unused_vertices(Fn,Vn)
-        constrained_segments = [[indFix1[c[1]]] for c in constrained_segments]  # Fix indices after point removal 
-
+
         # Remove 3 and 4 connected points
         E_uni = meshedges(Fn; unique_only=false)
         con_V2V = con_vertex_vertex(E_uni,Vn)
-        nCon = map(length,con_V2V)
-        Eb = boundaryedges(Fn)
-        indB = unique(reduce(vcat,Eb))
-
-        indLowCon = findall(nCon.<5)
-        indRemove = indLowCon[ [!in(i,indB) for i in indLowCon] ]
-
-        logicKeep = fill(true,length(Vn))
-        logicKeep[indRemove] .= false
-        indKeep = findall(logicKeep)
-        indFix = zeros(Int,length(Vn))
-        indFix[indKeep] = 1:length(indKeep)
-
-        Vn = Vn[indKeep] # Remove points 
-        constrained_segments = [[indFix[c[1]]] for c in constrained_segments] # Fix indices after point removal 
-
-        # Redo triangulation after points have been removed
-        TRn = triangulate(Vn; boundary_nodes=constrained_segments,delete_ghosts=true)
-        Fn = [TriangleFace{Int}(tr) for tr in each_solid_triangle(TRn)] 
-        Vn = get_points(TRn)
-
-        Fn,Vn,indFix = remove_unused_vertices(Fn,Vn)
+        nCon = map(length,con_V2V)        
+        indLowCon = findall(nCon.>0 .&& nCon.<5)
+        indConstrained = reduce(vcat,reduce(vcat,constrained_segments))
+        indRemove = setdiff(indLowCon,indConstrained)     
+        if !isempty(indRemove)
+            Vn,indFix = removepoints(Vn,indRemove)
+            constrained_segments = [[indFix[c[1]]] for c in constrained_segments]
+
+            # Redo triangulation after points have been removed
+            constrained_segments_ori = deepcopy(constrained_segments) # Clone since triangulate can add new constraint points
+            TRn = triangulate(Vn; boundary_nodes=constrained_segments,delete_ghosts=true)
+            Fn = [TriangleFace{Int}(tr) for tr in each_solid_triangle(TRn)] 
+            Vn = get_points(TRn)
+            constrained_segments = constrained_segments_ori
+        end    
+        Fn,Vn,indFix = remove_unused_vertices(Fn,Vn)    
 
         # Smoothen mesh using Laplacian smoothing
         Eb = boundaryedges(Fn)
@@ -4502,14 +4494,15 @@ function regiontrimesh(VT,R,P)
         end
         append!(V,Vn) # Append vertices 
         append!(F,Fn) # Append faces 
-        append!(C,fill(q,length(Fn))) # Append color data 
+        append!(C,fill(q,length(Fn))) # Append color data         
     end
     V,_,indMap = mergevertices(V; pointSpacing = mean(P))
     indexmap!(F,indMap)    
 
     return F,V,C 
 end
 
+
 """
     scalesimplex(F,V,s)
 Scales faces (or general simplices) wrt their centre. 
@@ -6509,4 +6502,13 @@ function perlin_noise(size_grid, sampleFactor; type=:Perlin)
         end
     end
     return M
+end
+
+function removepoints(V,indRemove)
+    n = length(V)        
+    indFix = zeros(Int,n)
+    indKeep = setdiff(1:n,indRemove)
+    indFix[indKeep] .= 1:length(indKeep)
+    V = V[indKeep]    
+    return V,indFix
 end
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -6722,3 +6722,16 @@ end
     @test_throws Exception perlin_noise(size_grid, sampleFactor; type=:wrong)
 end
 
+@testset "removepoints" verbose = true begin            
+    n = 15
+    V = Vector{Point{3,Float64}}(undef,n)
+    V_ori = deepcopy(V)
+    indRemove = [1,5,10] # Indices to remove
+    V,indFix = removepoints(V,indRemove)
+    indMap = [2,6,11,10]
+    indMapped = [1,4,8,0]    
+    @test V == V_ori[setdiff(collect(1:n),indRemove)]    
+    @test isa(V,Vector{Point{3,Float64}})    
+    @test length(V) == n-length(indRemove)
+    @test indFix[indMap] == indMapped
+end