Updated hemisphere for quads

examples/demo_hemisphere.jl

@@ -0,0 +1,41 @@
+using Comodo
+using GLMakie
+using GeometryBasics
+using Rotations
+using Statistics
+using LinearAlgebra
+
+
+# Demo for building a  hemissphere from a isosaheddron .
+
+# Example geometry for a sphere that is cut so some edges are boundary edges
+n = 1# Number of refinement steps of the geodesic sphere
+r = 1.0 # Sphere radius
+
+r = 1.0
+F1,V1 = hemisphere(1,r; face_type=:tri)
+F2,V2 = hemisphere(2,r; face_type=:tri)
+F3,V3 = hemisphere(3,r; face_type=:tri)
+
+F4,V4 = hemisphere(1,r; face_type=:quad)
+F5,V5 = hemisphere(2,r; face_type=:quad)
+F6,V6 = hemisphere(3,r; face_type=:quad)
+
+# Visualization
+fig = Figure(size=(1200,800))
+
+ax1 = Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=0")
+hp1 = poly!(ax1,GeometryBasics.Mesh(V1,F1), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+ax2 = Axis3(fig[1, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=2")
+hp2 = poly!(ax2,GeometryBasics.Mesh(V2,F2), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+ax3 = Axis3(fig[1, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:tri, n=3")
+hp4 = poly!(ax3,GeometryBasics.Mesh(V3,F3), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+
+ax4 = Axis3(fig[2, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=0")
+hp4 = poly!(ax4,GeometryBasics.Mesh(V4,F4), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+ax5 = Axis3(fig[2, 2], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=2")
+hp5 = poly!(ax5,GeometryBasics.Mesh(V5,F5), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+ax6 = Axis3(fig[2, 3], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "face_type=:quad, n=3")
+hp6 = poly!(ax6,GeometryBasics.Mesh(V6,F6), strokewidth=1,color=:white, strokecolor=:black, shading = FastShading, transparency=false)
+
+fig

src/functions.jl

@@ -1848,7 +1848,7 @@ end
 
 
 function pushtoradius_(v,r=1.0)    
-    return v.*(r/sqrt(sum(v.^2))) # Push to sphere       
+    return v.*(r/norm(v)) # Push to sphere       
 end
 
 
@@ -1894,47 +1894,46 @@ The algorithm starts with a regular icosahedron that is adjusted to generate a h
 Next this icosahedron is refined  `n` times, while nodes are pushed to a sphere surface with radius `r` at each
 iteration. 
 """
-function hemisphere(n::Int,r::T) where T <: Real
-    if n == 0 #creates a half a isosaheddron
-        F,V = geosphere(0,r) # Creating the faces and vertices of a full sphere
+function hemisphere(n::Int,r::T; face_type=:tri) where T <: Real    
+    if n<0
+        throw(ArgumentError("n should be >= 0"))
+    end
 
-        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
+    # Creating the faces and vertices of a full sphere mesh for "n=0"
+    if face_type == :tri        
+        F,V = geosphere(1,r) # A once refined icosahedron has an "equator" line
 
-        for (i,v) in enumerate(V)
-            if v[3]<0
-                v = Point3([v[1],v[2],0.0]) # Make z zero
-                V[i] = pushtoradius_(v,r) # Overwrite point        
-            end
-        end
-    elseif n>0
-        F,V = geosphere(1,r) # Creating the faces and vertices of a full sphere
-
-        ϕ = Base.MathConstants.golden # (1.0+sqrt(5.0))/2.0, Golden ratio
-        s = r/sqrt(ϕ + 2.0)
-        t = ϕ*s
-        α = atan(t/s) # angle needed for roating the sphere to 90 deg
-        Q = RotXYZ(0.0,α,0.0)
-        V = [GeometryBasics.Point{3, Float64}(Q*v) for v in V] 
-
-        VC = simplexcenter(F,V) # Finding triangle centre coordiantes
-        F = [F[i] for i in findall(map(v-> v[3]>=0,VC))] # Remove some faces at central coordinates z=0
-        F,V = remove_unused_vertices(F,V) # Cleanup/remove unused vertices after faces were removed
-
-        if n>1
-            for _ = 1:n-1
-                numPointsNow = length(V)
+        # Rotating sphere so a vertex "points north" and equator is in xy-plane            
+        Q = RotXYZ(0.0,atan(Base.MathConstants.golden),0.0) # Rotation matrix
+        V = [Point{3, Float64}(Q*v) for v in V] # Rotate vertices        
+    elseif face_type == :quad        
+        F,V = quadsphere(1,r) # A once refined cube has an "equator" line        
+    else
+        throw(ArgumentError("face_type should be :tri or :quad"))
+    end
+
+    # Now cut off bottom hemisphere    
+    searchTol = r./1000.0 # Tolerance for cropping hemisphere (safe, somewhat arbitrary if much smaller then mesh edge lengths)
+    F = [F[i] for i in findall(map(f-> mean([V[j][3] for j in f])>=-searchTol,F))] # Remove faces below equator
+    F,V = remove_unused_vertices(F,V) # Cleanup/remove unused vertices after faces were removed
+
+    # Refine sphere if needed 
+    if n>0 # If larger then refining is needed       
+        for _ = 1:n # Refine once n times
+            numPointsNow = length(V) # Number of points prior to refining
+            # Change refine behaviour based on face_type
+            if face_type == :tri        
                 F,V = subtri(F,V,1)
-                @inbounds for i in numPointsNow+1:length(V)
-                    V[i] = pushtoradius_(V[i],r) # Overwrite points
-                end
+            elseif face_type == :quad        
+                F,V = subquad(F,V,1)
             end
-        end
-
-    else
-        # error message here
+            # Now push newly introduced points to the sphere
+            @inbounds for i in numPointsNow+1:length(V)
+                V[i] = pushtoradius_(V[i],r) # Overwrite points
+            end
+        end        
     end
+
     return F,V
 end
 
@@ -2574,8 +2573,11 @@ function quadsphere(n::Int,r::T) where T <: Real
     V = coordinates(M)
     if n > 0
         for _ in 1:n
+            numPointsNow = length(V)
             F,V = subquad(F,V,1)
-            V = r .* (V ./ norm.(V))
+            @inbounds for i in numPointsNow+1:length(V)
+                V[i] = pushtoradius_(V[i],r) # Overwrite points
+            end
         end
     end
     return F,V