Simplify and correct the extrude function. (#1565)

* Add test

* Add getface and isclockwise funs.

* Check member variable exists before accessing it.

* sort_visible_faces() was really only working for triangles.

* Simplify and correct the extrude function.

src/psxy.jl

@@ -1465,14 +1465,16 @@ function sort_visible_faces(FV::GMTfv, azim, elev; del::Bool=true)::Tuple{GMTfv,
 	projs = Float64[]
 
 	!FV.bfculling && (del = false)		# Do not delete if bfculling is set to false (for example if FV is not closed)
+	first_face_vis = true
 	for k = 1:numel(FV.faces)			# Loop over number of face groups (we can have triangles, quads, etc)
-		n_faces, n_verts = size(FV.faces[k], 1), size(FV.verts, 2)	# Number of faces (polygons) and vertices of the polygons
-		tmp = zeros(n_verts, 3)
+		n_faces = size(FV.faces[k], 1)	# Number of faces (polygons)
+		this_face_nverts = size(FV.faces[k], 2)
+		tmp = zeros(this_face_nverts, 3)
 		del && (isVisible = fill(false, n_faces))
 		dists = NTuple{2,Float64}[]
 		_projs = Float64[]
 		for face = 1:n_faces
-			for c = 1:3, v = 1:n_verts						# Build the polygon from the FV collection
+			for c = 1:3, v = 1:this_face_nverts						# Build the polygon from the FV collection
 				tmp[v,c] = FV.verts[FV.faces[k][face,v], c]
 			end
 			this_proj = dot(facenorm(tmp, zfact=FV.zscale), view_vec)
@@ -1483,10 +1485,12 @@ function sort_visible_faces(FV::GMTfv, azim, elev; del::Bool=true)::Tuple{GMTfv,
 			end
 		end
 		data = del ? FV.faces[k][isVisible, :] : FV.faces[k]
+		isempty(data) && continue
 		ind  = sortperm(dists)
 		data = data[ind, :]
-		(k == 1) ? (FV.faces_view = [data]) : append!(FV.faces_view, [data])
-		projs = (k == 1) ? _projs[ind] : append!(projs, _projs[ind])
+		(first_face_vis) ? (FV.faces_view = [data]) : append!(FV.faces_view, [data])
+		projs = (first_face_vis) ? _projs[ind] : append!(projs, _projs[ind])
+		first_face_vis = false
 	end
 	sum(size.(FV.faces_view, 1)) < sum(size.(FV.faces, 1) / 3) &&
 		@warn("More than 2/3 of the faces found invisible. This often indicates that the Z and X,Y units are not the same. Consider using the `zscale` field of the `FV` input.")

src/solids.jl

@@ -339,7 +339,7 @@ end
 Create an extruded 2D/3D shape.
 
 ### Args
-- `shape`: The shape to extrude.
+- `shape`: The shape to extrude. It can be a 2D polygon or a 3D polygon defined by a Mx2 or Mx3 matrix or a `GMTdataset`
 - `h`: The height of the extrusion. Same units as in `shape`.
 
 ### Kwargs
@@ -354,25 +354,16 @@ function extrude(shape::Matrix{<:AbstractFloat}, h; base=0.0, closed=true)::GMTf
 		V = [shape; shape .+ [0 0 convert(eltype(shape), h)]]
 	end
 
-	if (h * facenorm(shape, normalize=false)[3] > 0)  v1, v2, f2, f4 = np:-1:1,     np+1:2np, 1, np
-	else                                              v1, v2, f2, f4 = 2np:-1:np+1, 1:np,     np, 1
+	if (h * facenorm(shape, normalize=false)[3] > 0)  v1, v2, f2, f4 = np:-1:1, np+1:2np,    1,  np		# Normal is pointing up
+	else                                              v1, v2, f2, f4 = 1:np,    2np:-1:np+1, np, 1		# Normal is pointing down
 	end
 
-	if (size(shape, 1) == 5 && shape[1,:] ≈ shape[end,:])		# When shape is a quad
-		F = (closed == 1) ? [[reshape([v1;], 1, :); reshape([v2;], 1, :); zeros(Int, np-1, 5)]] : [zeros(Int, np-1, 5)]
-		k = (closed == 1) ? 2 : 0		# Starting index of the faces vector that will contain the vertical faces
-		for n = 1:np-1					# Create vertical faces
-			k += 1 
-			F[1][k, 1], F[1][k, 2], F[1][k, 3], F[1][k, 4], F[1][k, 5] = n, n+f2, n+np+1, n+f4, n
-		end
-	else
-		if (closed == 1)  F = [[reshape([v1;], 1, :); reshape([v2;], 1, :)], zeros(Int, np-1, 5)]
-		else              F = [zeros(Int, np-1, 5)]
-		end
-		iF = (closed == 1) ? 2 : 1		# Index of the faces vector that will contain the vertical faces
-		for n = 1:np-1					# Create vertical faces
-			F[iF][n, 1], F[iF][n, 2], F[iF][n, 3], F[iF][n, 4], F[iF][n, 5] = n, n+f2, n+np+1, n+f4, n
-		end
+	if (closed == 1)  F = [reshape([v1;], 1, :), zeros(Int, np-1, 5), reshape([v2;], 1, :)]
+	else              F = [zeros(Int, np-1, 5)]
+	end
+	iF = (closed == 1) ? 2 : 1		# Index of the faces vector that will contain the vertical faces
+	for n = 1:np-1					# Create vertical faces
+		F[iF][n, 1], F[iF][n, 2], F[iF][n, 3], F[iF][n, 4], F[iF][n, 5] = n, n+f2, n+np+1, n+f4, n
 	end
 	fv2fv(F, V; bfculling=(closed == 1))
 end
@@ -409,7 +400,7 @@ Return a Faces-Vertices dataset.
 function cylinder(r, h; base=0.0, center=(0.0, 0.0, 0.0), closed=true, geog=false, unit="m", np=36)::GMTfv
 	h0 = (base != 0.0) ? base : length(center) == 3 ? center[3] : 0.0
 	if (geog == 1)
-		xy = circgeo(center[1], center[2]; radius=r, unit=unit)
+		xy::Matrix{Float64} = circgeo(center[1], center[2]; radius=r, unit=unit)
 	else
 		t = linspace(0, 2pi, np)
 		xy = [(center[1] .+ r * cos.(t)) (center[2] .+ r * sin.(t))]

src/utils.jl

@@ -990,6 +990,38 @@ Return the number of rows, columns and segments in the dataset `D`.
 """
 getsize(D::GDtype) = isa(D, GMTdataset) ? (size(D.data,1), size(D.data,2), 1) : (size(D[1].data,1), size(D[1].data,2), size(D,3))
 
+# ------------------------------------------------------------------------------------------------------
+"""
+	getface(FV::GMTfv, face=1, n=1; view=false) -> Matrix{Float64}
+
+Return the n-th face of the group number `face`. If `view` is true, return the viewable face.
+
+### Example
+```julia
+FV = cylinder(1.0, 4.0, np=5)
+getface(FV, 1, 1, view=true)
+```
+"""
+function getface(FV::GMTfv, face=1, n=1; view=false)
+	(view == 1) ? FV.verts[FV.faces_view[face][n,:],:] : FV.verts[FV.faces[face][n,:],:]
+end
+
+# ------------------------------------------------------------------------------------------------------
+"""
+    isclockwise(poly::Matrix{<:AbstractFloat}, view=(0.0,0.0,1.0)) -> Bool
+
+Return true if the 2D/3D `poly` is clockwise when seen from the `view` direction.
+
+### Example
+```julia
+poly = [0.0 0.0 0.0; 0.0 0.0 1.0; 0.0 1.0 1.0; 0.0 1.0 0.0; 0.0 0.0 0.0]
+isclockwise(poly, (1.0,0.1,0.0))
+```
+"""
+function isclockwise(poly::Matrix{<:AbstractFloat}, view=(0.0,0.0,1.0))
+	dot(facenorm(poly, normalize=false), [view[1], view[2], view[3]]) <= 0.0
+end
+
 # ------------------------------------------------------------------------------------------------------
 """
     settimecol!(D::GDtype, Tcol)

src/utils_types.jl

@@ -936,7 +936,7 @@ function rasters2grid(arg; scale::Real=1f0, offset::Real=0f0)::GMTgrid
 			(for k = 1:numel(names) names[k] = names[k][1:10] end)
 	end
 
-	dic = arg.metadata.val
+	dic = !isempty(arg.metadata) ? arg.metadata.val : Dict()
 	(scale  == 1) && (scale  = get(dic, "scale", 1.0f0))
 	(offset == 0) && (offset = get(dic, "offset", 0.0f0))
 	z_unit  = get(dic, "units", "")

test/test_solidos.jl

@@ -13,4 +13,6 @@
 	D  = GMT.replicant(FV, replicate=rand(5,3)*100);
 
 	FV = gmtread("file.obj");
+	FV = cylinder(1,4, np=5);
+
 end